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SIMULATION OF THE FOKKER-PIANCK EQUATION BY RANDOM WALKS 
OF TEST PARTICLES IN VELOCITY SPACE WITH APPLICATION 

TO MAGNETIC MIRROR SYSTEMS 
by Gerald W. Englert 
Lewis Research Center 

SUMMARY 

The structure of an analytical expression obtained from a random-walk process was 
compared with that of the Fokker-Planck equation. The step sizes and probabilities of 
taking steps in various directions were related to the coefficients of the Fokker-Planck 
eqtiation. i^herical coordinates with azimuthal symmetry about the polar axis were 
used. 

The technique was applied to a magnetic-mirror system for confining charged par¬ 
ticles. The cases selected have approximate analjrtical solution via the Fokker-Planck 
equation. Velocity distributions inside and end losses from this system were determined 
by the random-walk procedure and compared with results of the analytical solutions. 

The effect of reducing the maximum impact parameter cutoff distance well below the 
Debye length was studied. This reduction increased step sizes and decreased collision 
frequency in such a manner that the probability of a test particle walking to specified 
locations in velocity space remained approximately constant. A sampling of random 
numbers from the lower impact-parameter range greatly reduced computing time and 
yet permitted reliable distributions within 10- to 20-percent accuracy. 

The random-walk results were, in turn, permitted an appraisal of some simplifying 
assumptions made in the analytical solution of the Fokker-Planck equation. An iteration 
procedure was used to find a self-consistent solution of the distribution of test particles 
walking throi^h an ensemble of field particles, on which the Fokker-Planck coefficients 
were based. Analytical solutions assume that the distribution of test particles and loss 
rates are insensitive to assumed field-particle distributions. For the cases studied this 
was found to be a good assumption. Two steps in the iteration process made the 
random-walk solutions self-consistent within the scattering of the data for a sampling 
size of 500 test particles. The popular assumption of separability of the distribution 
function was also verified. 

The random-walk method should apply with little additional difficulty to a number of 
problems that cannot be solved analytically. 



INTRODUCTION 


The Fol^er-Planck equation, in general, describes the time development of a 
Markov process. Such a process is characteristic of the nature of classical collisions 
where each event depends on the present conditions and is independent of the past 
(ref. 1, p. 157, and ref. 2, p. 369). Written in velocity space, this equation provides 
a popular approximation to the collisional terms in the Boltzmann equation (refs. 3 to 5). 
Included in this approximation are first and second moments of velocity increments de¬ 
scribing two-body encoxmters (p. 75 of ref. 6). These moments enter as coefficients 
of the Fokker-Planck equation and permit it to describe changes in the velocity distri¬ 
bution function resulting from influence of dynamic friction and dispersion (refs. 

7 and 8). In general, the velocity distribution function is used to weight the moments 
of velocity increments, making the Fokker-Planck equation nonlinear and very difficult 
to solve for many problems of interest. 

Another approach to collision problems is through the study of random walks of 
test particles. This approach replaces the mathematical complexity with a relatively 
simple but repetitious process, which is greatly facilitated by high-speed computers. 
Usually the walks of a large sample of test particles through an ensemble of field par¬ 
ticles must be studied to determine how they distribute in velocity space. For a self- 
consistent solution, the resultii^ test-particle distribution should coincide with the 
field-particle distribution, requirir^ an iteration process. 

A random walk will be defined herein as a process in which step sizes and proba¬ 
bilities of taking various steps are average values (depending on field-particle distribu¬ 
tions and test-particle location in velocity space for the case at hand). In contrast to 
this procedure is the more popular Monte Carlo method. This technique makes random 
selections from distributions of pertinent variables (such as impact parameter and rela¬ 
tive velocity) at each collision and performs an interaction calculation at each step to 
determine the resulting test-particle location (refs. 9 and 10). Such a detailed calcula¬ 
tion at each step is excellent for physical and mathematical clarity, but would involve 
an excessive amount of computer time for cases involvii^ a very large number of en¬ 
counters. On the other hand some physical clarity can be easily lost in the representa¬ 
tion by random walks, especially in determining the probability of takir^ steps in vari¬ 
ous directions. The ability to relate step sizes and probabilities to the well explored 
coefficients of the Fokker-Planck equation would thus be very useful. 

The random-walk and the Fokker-Planck concepts depend primarily on the same 
combinatory laws of probability; however, the random walk as depicted herein is re¬ 
stricted to steps on a grid. A detailed generation of the Fokker-Planck equation from a 
random-walk procedure is available in the literature only for the one-dimensional prob¬ 
lem with constant coefficients. Reference 2 (ch. 14) shows that for this case the proper 
step size of a random walk is dependent on the second moment and that the probabilities 
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of taking steps in positive and negative directions depend on both first and second mo¬ 
ments. This derivation is extended herein to problems involving three dimensions and 
variable coefficients. In this more general case the step sizes and probabilities are 
again determined in terms of Fokker-Planck coefficients. The Fokker-Planck coeffi¬ 
cients are evaluated for the binary Coulomb encounters used in plasma physics. 

A main difficulty in obtainii^ final results is due to the loi^ range nature of Coulomb 
encounters, which makes the step sizes of the walks extremely small. For example, at 

7 

typical experimental or laboratory conditions, a particle may average 10 encoxinters 

(steps) before it reaches a 90° deflection from its initial direction. At thermonuclear 

1 fi 

conditions this number easily reaches 10 steps. Tracing such a large number of 
steps for a representative number of test particles would require an excessive amount 
of time even for today's fastest computers. A selective sampling technique of some 
nature must be therefore employed to the random-walk procedure. 

The coefficients of the Fokker-Planck equation are insensitive to the cutoff assumed 
for the collisional impact parameter (refs. 11 and 12). As the maximum impact param¬ 
eter is reduced, the average step size increases, but the number of encounters (steps) 
per unit time decreases. With a low impact-parameter cutoff distance a test particle 
would take a relatively small number of large steps. On the average, however, it would 
reach a specified boundary in velocity space in about the same time as in a walk based 
on a higher impact-parameter cutoff and requiring a large number of small steps. The 
computii^ time, however, varies inversely as about the square of the step size and is 
thus much less for walks based on lower impact-parameter cutoff distances. 

To study results of this sampling technique, it was applied to the end-loss problem 
of magnetic-mirror systems. The scattering of charged particles into a loss cone has 
received considerable attention over a period of years. However, because of the com¬ 
plexity of the Fokker-Planck equation, the rate of such loss has been determined ana¬ 
lytically for only the simplest of cases (refs. 5 and 3). Numerical solutions (ref. 13) 
are available for only a restrictive set of initial and boundary conditions. Desired also 
are particle distributions inside the mirror system for use in stability studies (ref. 14). 

In summarizing, the attempt herein is to exchange a difficult analytical problem 
with a fairly straightforward computing procedure. Computing time is made reasonable 
by a sampling technique. Considerable physical detail can be incorporated into such a 
procedure with little increase of mathematical complexity. 


ANALYSIS 

The general form of a random-walk equation will be derived for walking on a co¬ 
ordinate grid in N dimensions. (One-dimensional walks with constant step sizes and 
probabilities are treated in refs. 2 and 8. Ref. 2 writes the results in the form of a 



Fokker-Planck differential equation. Ref. 8 uses the random walk to derive solutions 
to some boundary value problems, but does not arrive at a differential equation as an 
intermediate step. The Fokker-Planck equation is also derived in ref. 8, but not from 
the concept of a random walk on a grid.) It is not necessary for the grid to be orthog- 
ornal or have constant spacing. Comparison will then be made with the Fokker-Planck 
equation applied to inverse square collisions in velocity space. The results will be 
adapted to a study of charged particles in a magnetic-mirror system. 

A list of symbols is given in appendix A. The International System of Units (SI) 
will be used throughout with the exception of temperature T being reported in keV and 
the corresponding Boltzmann constant k in joules per keV. 


Development of the Random-Walk Equation 

Let 1 2 , lo’ ■ * ■ ’ In independent coordinates that span the space of interest. 
Let Pj and q^ be the probabilities that the i*" coordinate will increase or decrease 
during the course of a step. To identify these conditional probabilities, only their loca¬ 
tion at the start of a step will be labeled in their arguments. For example, compared 
with the more conventional nomenclature. 


^ 2 ’ * * •’ 


means 


and 




^ 2 ’ • • • > 


means 


Ql(^l ” ^2^ • • •> ^2’ • * 

Let the probability of a step being in the ii*^” direction be 1/N; that is, 

Pi(ll> ^2’ * * •’ ^2’ ' • —— (1) 

N 

Let l2’ • • •» denote the chance of a particle being at location 
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(a) Steps along component direc¬ 
tions; number of components 
changed per encounter, 1. 



(b) Steps along diagonals; number of com¬ 
ponents changed per encounter, 2 . 


Figure 1. - Examples of random walks in two- 
dimensional Cartesian coordinates. 


J, ^ 2 ’ •• •’ interval between 

starts of the n and n + t) steps. 

Consider the case for t) = 1, illustrated in figure 1(a). The change in the number of 
particles at (|j, ^ 2 ’ • • • > during At is the probability of reaching (4j, Ig* • • •> 

from the closest neighbors minus the chance of a particle leaving (^j, • • •» 

and going to a nearest neighbor. Since a step of finite size is taken each time increment, 
there is no chance of a particle standing still. The conservation of particles at (^j, 



can thus be written 
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^2’ * * ’ ' ^2’ * * ’ ’ ^2' * ’ *' ^2' ‘ * *' 

+ Qlfel+^^1> ^2’ * * ‘J 1 ^2’ • • •> 

+ P2(^iJ ^2 ” ^^2’ * * *» ^2 ~ ^^2’ * • •’ 

+ ^2^^!’ ^2 ^ ^^2’ * * ’’ ^2 ^ "^^2^ * • •> ^jq^) 

+ * • • ?2’ * ' ‘^ % " ?2’ * * •' ^N ’ 

+ ^2’ • * *’ ^N ^2’ ’ * ' ^ ^N 

- Pl(^lj ^2’ ■ ■ ‘ ’ ^2' * * * ’ ^2’ ' ■ *’ ^2’ * * *> 

- P2(?l^ ^2" * ■ *' ^ 2 ’ * * ■" ~ ^ 2 ’ * ' ^2' * * *’ 

- • • • ^ P-^i^ii ^ 2 ' * * *' ^ 2 ^ ' * ■ ’ ^ 2 ’ * * * ’ ^ 2 ’ • ’ *» %) 

where use is made of the first and second laws of composition of probability. 
Multiplying through by 2N and adding and subtracting u^^ terms results in 

• * • > ^2’ ' * •’ %)] = 42> ■ • •’ ^2’ * ' ' ’ 

- [1 - NpjClj - A|j, ^2j • • - j ~ ^^1^ ^2^ ’ ‘ •r ^■^) 

■** Nqj(|^ + A|^, ^25 • • • > ^ ^2’ • • • » 

- [1 - Nqj + A| j; Ig; • ' - ^ ^2^ * • w Ijq) 

^2’ ■ ■ • ’ ^2’ ■ * *’ ^^1’ ^2’ • • •’ 

+ Np2(^l, ^2 " ^^2’ • * *’ ^2 “ ^^2^ ' * ’ ’ 

-[1-Np2(^j, ^2 ■ -^^2’ • • -^ ^2 ■ ^^2’ ' ’ ’ ’ 

+ Nq2Ui, ^2‘*'^^2’ * • •» ^2’ ' ‘ ’ ’ 

-[1 - Nq2(l2, ^2"^ ^^2' ' ‘f ^2"^ ^^2^ • • ., Ij^) 

+ ^2 " ^^2' • • •» ^2^ * ■ ‘ ’ ^2 ^^2’ • ' * > 

+ . . , + Np^dj, ^2? * • • > ^N ” ^2’ ■ ■ ■’ *“ 

- [1 - Npj^(|j, ^2 j • - • > ~ ^2’ ‘ ■ ■ ’ ^N “ 

+ Nqj^dj^, ^2’ * ‘ * ’ ^N ^2’ * ' ' ’ ^N 

- [1 - Nqj^(^j, ^2y • * • > 4jij + ^2’ * * * ’ ^ 

+ ^2’ * • •» ^N ■ ^2’ ’ ’ ’^ ^2» * • *’ ^N ■*' 
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Using equation (1), combining and terms of like arguments and i subscripts, 
collecting u^ terms and finally writing the result in finite difference notation yields 


2^An = 2^At=-2{6i[(Pi-qi)uJ+62[(P2-q2)u2]+- • • + j} + ^ ' • • + ®NnV 

An At '• JN ^ 

which approaches the limit 


N 

3u_ V / ^BPj - qj)^ ^^i , 1 e^u 

9t a|. At 2g^2NAt J 

i=l ^ 

as step size (grid spacing) becomes arbitrarily small. 

This is the more conventional approach to a random-walk expression. However, 
next consider the case where steps are taken in groups of tj = N; one along each coor¬ 
dinate direction. (This corresponds to the later physical description of a test particle 
in polar coordinates where its magnitude and direction (v, 9, and cp) change each en¬ 
counter. ) This corresponds to steps across the diagonals of N-dimensional lattices 
which will herein be called path lines. This is illustrated for the case of N = 2 on 
figure 1 (b). 

Repeating this procedure with 77 = N gives 

9 [ n(p^ - q^u] ^ ^ 1 a^u 

a|. At 2 ..2 At 

1 9 |. 

This same result could be obtained by letting 

Pj(^j> • • • > ^2’ ■ • ■ ’ i — 1, 2, ...,N 

in the derivation regardless of N, and letting ?? = 1; that is, there is a sure chance of 
stepping either in the plus or minus direction of each of the components for each time 
increment At. 

Equation (2) is in the form of the Fokker-Planck equation except that no mixed 
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partial derivatives are present. Generation of these terms does not appear possible 
when the walks are restricted to the nearest neighbors on a grid. Their magnitude will 
be discussed later. 

Note that p. - q. is like a bias in the i direction and that the terms involving 
first-order derivatives could be, for example, related to friction. To keep the process 
from being degenerate, p^ - q. must approach zero as step size approaches zero (p. 324 
of ref. 2). The terms involving second derivatives are similar to diffusion expressions 
with (A|.)^/At similar to a diffusion coefficient. 

When a stochastic process is beir^ simulated by a random walk on a grid, the step 
sizes are often statistical averages. It can be deduced from references 2, 8, 15, or 16 
that second moments determine the proper grid spacing. Consider, for example, the 
application discussed in the next section. In any small region of velocity space the test 
particle taking a random walk encounters many other field or background particles, one 
at a time. A field particle can be in any accessible region of velocity space as it starts 
its encounter with a test particle (its probable location is defined by appropriate distri¬ 
bution functions). The grid spacings are thus based on statistical averages of encounters 
with field particles (fig. 2), these averages nevertheless being dependent on test- 
particle location. 


Pi and qj 



Probabilities of increase and decrease, 
respectively, of i^^ component 
(i = 1, 2, 3) 

i^*^ component of step size 


^3 
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Fokker-Planck Equation 


The widely accepted Fokker-Planck equation for studying plasmas in N dimen¬ 
sional Cartesian coordinates is 



— (f<Av.)) +-- - - 

3v. ^ 2 0v. 9v. 

^ J 


(f<Av. AVj)) 


where (9f/0t) denotes the change of probability density (distribution function of test 
particles) with respect to time due to collisions. The coefficients (Av.) and (Av. Av.) 

J 

are averages over distributions of field particles and scattering angle (ref. 11). The 
angular brackets signify a time average; that is, (Av^) is the average increment per 
unit time of the i’^^ component of velocity. The average value per collision, denoted 
by a bar, is equal to the time average divided by the collision frequency p>‘, for example, 
Av^ = (Av^}/u. 

The Fokker-Planck equation is transformed to spherical coordinates with symmetry 
about the polar axis in appendix B. This is essentially the same coordinate system as 
used in parts IV and V of reference 11. In reference 11, cos 9 is used in place of the 
polar angle 0 used herein. The assumption of aximuthal symmetry is not too restric¬ 
tive for many problems in plasma physics and serves mainly to reduce cumbersome ex¬ 
pressions. The analysis follows in a straightforward manner without it. In this coor¬ 
dinate system 


(• 


= . i. A (v2f( AV» - ^ A (f sin 0< A0» + J_ ^ sin 9((Ae)^ 


Bt/c v 2 3v 


sin 0 39 


2v^ 3v^ ^ ^ 


30' 


^ 1 a^v^f sin 9<Av Ae) _ J._ _a_ |v%[( (A0)2) + Sin2e( (A(p)2)]| 

v2 sin 0 30 2 v 2 3^ 


2 v sin 


I - —if sin 0r2(A0 Av) - V sin 0 cos 0( (A<;p)^)])[ 

Bin 0 90 1 ^ 



Using the identities 


8Vf<(Av)^) ^ ^ 2 ± (vh 9<(Av)^)\ _ ^2^ 9^((Av)^) 

av2 av2 9v V av / 3^2 


and 


ah sin eHAe )h , ah Bin a 2 A 4 si„ « i<(M!>'\ . (sm e 

dO^ de^ ' 35 / gg2 


the Fokker-Planck equation can be written as 


(• 


2//a—\2\ ^2 




K n 3^2 


=-A±U <AV> -i<(^+I<(A0)2) +Zsin2fl((A,^)2)1 

gg2 J ^2 dvl av 2 2 J 

--J-Aff Sin0/<A5) - 9 . COS ^ ((a^)2))] 

sin 0 30 L \ V 2 / J 2 • 


3^v^f sin 0 ( Av A0) 
3v 30 


<(Av) 2> a2fv2 ^ ({^9)^) dH sin 9 
2 v 2 3 v 2 2 sine g^2 


2\ ^2r 


(3) 


Correspondence Between the Random-Walk and Fokker-Planck Equations 

The random-walk expression (eq. (2)) describes the average results of a large num¬ 
ber of walks, or a large number of walks gives a numerical solution to the random-walk 
equation. An attempt will now be made to write equation (2) in a form in close agree¬ 
ment with equation (3) to see whether random walks can provide solutions to the Fokker- 
Planck equation. 

In equation (2) let N = 2, = v, At = 1/v (where 

i> = i/(v, 0) = collision frequency). By equation (1) Pj + q^ = 1/2 and P 2 + q 2 = 1/2* Let 
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2(Pj^ - q^) = 


1 _9((Av)!) ^ V ^ V sin^ ^^^^2 
1 / 0v 2 2 - 

V 


and 


2(p 2 - qg) = 


A9 - - 

_ - V 


1 9<(A0)^> 


A0 Av , sin 9 cos 9 /. „\2 
-+- (Acp) 

V 2 J 


(A0) 


All = V (Av) 


A^^ = \{A9) 


u = fv^sin 9 


The random-walk equation then becomes 


(M.) -1 (Av> - 3((^1 ^Z((a9)2) +1LEH1:£<(a^)2>’ 
Vat/j. 2 [V 3v 2 2 J 

+ f<Ae> - c° .g.£((A<p)^)) /- 

\ 30 V 2 / Y 


1 3((Av)^) ^ 1 dy 

((Av)2> 


1 3((A0)^) , 1 3v 

((A0)2> 9® '' 


= -i.AUf(Av> - 3<(^+Z<(A0)2> ^vgLni9 <(A<,)2)\1 
^2 0vL \ 0v 2 2 /J 

-_i_±rf sin0f<A0) - 3((A0)^) <A9 Av) , sin 0 cosj^ <(A.v)!> jV|+i 

sinSaaL \ 90 V 2 /J 2 2 , 2 2 sin 0 gg2 
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which differs from the Fokker-Planck equation by the terms involved with f/2 on the 
left side of the equal sign, and by the absence of the mixed derivative term on the right 
side. 

To study these equations further and to perform random walks, one must specify 
the initial and/or boundary conditions. To evaluate the Fokker-Planck coefficients, one 
must specify the type of encoimter. The differences between equations (3) and (5) will be 
appraised for specific conditions. 


Physical Application 

The preceding results will be applied to magnetic mirror systems used to confine 
ensembles of electrically charged particles. The Fokker-Planck coefficients will thus 
be based on Coulomb encounters. The random walks will determine the ion distribu¬ 
tions inside and the end loss from a magnetic mirror system. Such a system is sketched 
in figure 3 and discussed, along with simplifying assumptions, in appendix C. The ve¬ 
locity space of prime interest can be best described in spherical coordinates (fig. 4). 

The polar axis is alined with the magnetic field about which there is axial symmetry. 

This B field is assumed to be uniform over the central portion of the physical space. 
The increased magnitude of B at the mirrors enters the problem by describing a loss- 
cone boundary condition in velocity space. The random walks terminate at this boundary 
as if it were an absorbing wall. These assumptions result in miiformity of f in <p and 
in coordinate space and, thus, reduce a problem in six-dimensional phase space to two 
dimensions in velocity space. 


Mirror magnet 



Figure 3. ~ Magnetic mirror system. 
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The usual simplifying assumptions applied to the analysis of magnetic mirror sys¬ 
tems reduce the Boltzmann equation to 


= -s (6) 

c 

for steady state (appendix C). The symbol s represents a sink (or source) density in 
velocity space. For a steady state the rate of particles being injected into the system 
equals the rate of particles entering the loss cone. Loss rates from the system can 
therefore be expressed as 



n = 




where integration is over the velocity space of the injection distribution. The injection 
distribution provides the initial condition for the random walk. 

To correct the random-walk results for the difference in the expressions on the 
left sides of equations (3) and (5), the following equation can be used: 


[ 
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„ _n f re^((Av)^) , 3^<(A9)^) 

^corrected “random walk / 9 1 9 9 


A 


(Av) - 9<iMl+Z<(A9)2) 

av 2 


+ vsin!0<(A<pn 


2vW 1 3 <(Av)2) ^ 1 du 

<(Av)2> 3^ 3'^> 


f(A9) _ 9((A9)^) _ (A9 Av) sin 9 cos 9 ((a..)4[_1 

\ 3® 2 /\<(A9)2) 


( 8 ) 


3<(A9)^) ^_1 V 
30 ly 30) 


►dv 


The magnitude of this integral was negligible for cases of most interest. (This will 
be discussed later.) When the mixed-derivative term of the Fokker-Planck equation is 
negligible, such as for near isotropic field distributions, "corrected 
determined from the Fokker-Planck equation. 


Fokker-Planck Coefficients 


Expressions for the Fokker-Planck coefficients for the preceding physical applica¬ 
tion will now be supplied so that the step sizes and probabilities of equations (4) can be 
determined. 

The general integral expressions for these coefficients are derived in appendix D 
and can be expressed as 



<A0) 


= . 2r /* ^ - cos 9 sin cos(^ - <p^)] 


V + - 2vv^ cos 


3/2 


f(Vb)dVb (9b) 


•r. I vv, (v - V, cos 3) [sin 9 cos 9. - cos 9 sin 0^ cos(^ - 

AV) 4 Jf ;. -3/, «v.)av. 



^ ^ (v^ '^b “ ^r) 


(9c) 
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<(Av)^> = r 



(v - cos S^) 


(v^ + vj - 2^rv^ cos 3^) ' (v^ + " 2w^j cos 5 J 


3/2 


f(v^,)dv^, 


(9d) 


{(AS)2)= X 

A 


v^v^jsin 6 cos 0^ - cos S sin 0^ cos((p - 


^ - 2w^ cos Oj.) ^ 


(■ 


2 2 

V + - 2w^ cos 


3/2 


•f(vj,)dv^j (9e) 


({Acp) ) = 


2x r 



2 vvjj sin cos {cp - (p^^) 


V - 


sin 9 


sin^0 f /, 2 2 o Q 

(v + Vs - 2w^ cos 5 


■) 


1/2 


f(Vb)dv^ 


(9f) 


4ri3^ / m 

z/ = ^ ' 


9 In ^3 \2kT^ 


j y* (v^ + -2w^ cos Oj.) ^ f(v^,)dv^ 


(9g) 


where ^ is the ratio of maximum to minimum impact parameter and 


r = /3 _ 0. 24 In ;3 

.1 2 2 .2 

47re m A 


is the interaction parameter. 

In the usual problems of interest the field distribution f(^) is not known. Also, a 
test particle is usually just a background particle with a special label. Tracing the 
paths of a large enoi^h sample should determine test-particle distributions, which for 
a self-consistent solution should be the same as the fie Id-particle distribution. Such a 
solution can be approached by an iteration procedure. A background distribution is first 
assumed; the coefficients are calculated; and the paths of a representative number of 
test particles are traced to determine their distribution. The distribution of test par¬ 
ticles in the step of this iteration are used for the field-particle distribution in the 
I + \ step. This process is continued until the field- and test-particle distributions 
agree within the desired accuracy. The accuracy, or amount of scatter of the final data, 
in the example worked herein, was mainly set by the number of test particles used. 

In the first step of the iteration process, use of a Maxwellian distribution is con- 
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sistent with the theory of references 5 and 17. It is argued in these references that loss 
rates should be insensitive to the assumed field-particle distributions. For purposes 
herein, this may be considered a first-order solution. The insensitivity to field distri¬ 
bution may indicate convergence of the iteration process in a small number of steps. 

Maxwellian dis tribution of field part i cles . - A Maxwellian field distribution should 
provide an especially good first approximation to plasmas which are in near equilibrium 
conditions. An example of such could be inside mirror systems having very high mirror 
ratios so that field defects due to loss cones would be small. 

Since a Maxwellian distribution is isotropic, all first moments of 9 are zero. (An 
isotropic distribution corresponds to the diffusion approximation of ref, 5 and to the use 
of effective Rosenbluth potentials in ref. 11.) Thus {A9 Av) is zero, and the mixed de¬ 
rivative term is absent from equation (3). Use of 


f(v^)dv^ = % 


m 


27rkT 


,3/2 2 - mv^/2kT, 
% ® 


sin 


d^^^ dv^ 


b- 


( 10 ) 


and integrating between the limits 

0 < Vi < °o 

~ b “ 

0 < < TT 

0 < S: 21T 

reduces the expressions for the Fokker-Planck coefficients to the following (see appen¬ 
dix D): 


<Av> 


2 n^r 


4 . 


m 


2jTkT, 


V e 


-mv^/2kT^ 



<(Av)2) = - 



mv 


(Av) 


(11a) 


(lib) 


< (a@)2) = _L. 


- » Trim 


2kT, j -mv^/2kTj^ 


jrm V 



(11c) 
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<(A<P)2> =^L_ <(A0)2> (lid) 

sin^0 



(lie) 

Using equations (11) in (3) finally reduces the Fokker-Planck equation to 


/M\ .i<i^ + v{(A0)2>)] . <(Ag)^> -9,(fcos 

\dt/^ 2 3^2 ^2avL V av /J 2 sin0 30 9^2 

+ a^f sin 0 

2 sin 0 9^2 (12) 


An interesting relation among the coefficients of equation (11) is 


J- ± fv2<(Av)2>) - v<(A9)2) ^ 1 ^3) 

2v 2 av V ^ 2 

Use of equations (6) and (13) reduces (12) to equation (18. 7) of reference 5. 

Equation (12) is especially useful because there is an approximate analytical solu¬ 
tion to it in reference 5. The random-walk results will later be compared numerically 
to this solution. 

Using equations (11) reduces the random-walk equation (5) to 


(K] ,i[/iAv) ■ j<(^^v<(Ag)2))/_ L _9<i^^i 

\dt/^ 2[\ av '\((av)2> 3v 

-.J_Arv2f/<Av) - 9<(^ + v<(Ae)2)^l (14) 

^2avL V 3V /J 2 Sins as av^ 2 Sins gg2 
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In terms of dimensionless variables the average value per collision of the variables 
in equation (11) can be conveniently written as 


AV = 


(Av) _ 9 hi ^ 


9) 



(15a) 


(AV)^ = ^ ^ AV 


2v 


2VE’ 


o 


(15b) 


2 9 In /3 


{Aer 


4 v^/3^ 


9) 



(15c) 


(Acp)^ = _J_ (A0)2 
sin^0 


(15d) 


3/2i^l/2 

—^-= 0. 518X10"^^ 


4 


1 + 


2V^E’ 



rf 


/ I TP * \ 

Ufs ) + 1 

\ 1 T^/ (nF' 
\ b/ ^ 


-V^E^/T^, 


(15e) 


where 




E = i mv^ 
0 2 0 


E’ = — 

o 


i3 » 1 


18 



and 
enced to 


is some reference test-particle velocity. The quantity E' is always refer- 


and velocity v is referenced to v ; thus, v is actually referenced to 


through E^. If E*/T^ = 3/2, then by the kinetic definition of temperature E^ 

(p. 37 of ref. 18) is the average field-particle energy, and V is the ratio of the test- 
particle velocity to the rms field-particle velocity. Thus for most applications 
herein E*/T|^ was set eqml to 3/2. When studyii^ the case of injection of a mono- 
energetic beam of test particles into a background ensemble, however, it was conven¬ 
ient to set Eq equal to the beam-particle energy. Thus for example, when the beam 
particles had an energy of 10 times the average field-particle energy, E /T^ was set 
equal to 15, making the initial value o f V eq ual to 1. _ 

Plots of generalized step sizes V (AV)^ ^jS^E^/T^y^n^) and vV 

|/ln /s) for a Maxwellian field distribution are shown in figure 5. These curves were ob¬ 
tained by use of equations (15a) to (15c). 

Computir^ tim e and impac t pa rameter cu toff distance . - Setting the maximum im¬ 
pact parameter at the Debye length and the minimum value at a distance to cause a 90 
deflection results in 


o 


^Debye 


0. 49x10^® Te^^ 


2„l/2 


(16) 


Z"n 
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-3 

for (in Kev) and n^^ (in m“ ). Numerical values are given in table 5. 1 of refer¬ 
ence 12. 

If the maximum impact parameter is cutoff at the Debye leng th as in e quation (16), 
the average deflection which a particle undergoes per collision is extremely 

small, as can be deduced from figure 5. Thus the corresponding computing time re¬ 
quired to trace representative samples of test particles would be prohibitively large. It 
was observed from the elementary one-dimensional (and nonbiased) random walks of 
reference 8 and also from preliminary computer results that the average number of 
steps to reach a distance m steps away from a starting point was nearly proportional 
to m . The time to take such a walk would thus be nearly proportional to m /p. For 
a random walk in the 9 direction the distance reached per unit time would correspond 
to u{A9)^, which is simply the diffusion coefficient ((A0)^) . 

Solutions to processes described by the Fokker-Planck equation are dependent on 
impact parameter cutoff distance only through the Fokker-Planck coefficients. As shown 
by equations (9) and (11), these coefficients are only weakly dependent on impact param¬ 
eter ratio It thus appears that a particle taking the fewer number of larger steps 
determined by a lower p cutoff would reach given locations in velocity space in about 
the same time as if it took the larger number of smaller steps determined by a higher (B 
cutoff point. The final walk time (loss rate) would possible be weighted by a function of 

As part of the following computer study, the effect of on loss rates from and 
distributions inside magnetic mirror systems was determined. Cases were selected 
that could be solved analytically for comparison with the random-walk results. 


Numerical Procedure 

The procedure by which the random walks were simulated by the computer and the 
manner in which velocity distributions inside and end losses from the magnetic mirror 
system were obtained from the walks will now be described. 

A computer flow diagram and a representative FORTRAN program are presented 
in appendix E. Either the Lewis IBM 7044-7094 or 7040-7094 direct coupled systems 
were used to perform the calculations. 

Initial conditions for each random walk were determined by selecting at random 
from prescribed injection distributions. The walks were terminated at the loss-cone 
boundaries 9 and tt - 9 where 
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( 17 ) 



and B^/B^ is the ratio of the magnetic field at the mirror to that in the central uniform 
region (ref. 12). 

In any of the steady-state models selected for study herein, the average number of 
particles injected per unit time equals the rate at which particles enter the loss cone. 
The walks were studied one at a time, and a large enough number of them were traced 
to form a representative statistical sample. 

The velocity space to which the test particles had access was divided into zones . 

F or each of these zones values of probabilities p^ and Pg, step sizes ^{{Av)^)/v and 
y({A9)^}/v, and collision frequency were computed for the field distribution of inter- 
est. These quantities were stored in matrix form and called for as the test particle 
progressed about the velocity field, as described in appendix E. 

During each step of each walk, a random number was selected for each of the two 
components (V and 9) and compared with probabilities Pj and pg, respectively, to 
determine step direction. These random numbers were selected from a set of numbers 
uniformly distributed over the interval from 0 to 1 (ref. 10). Whether the V compo¬ 
nent of the step were in the positive or negative direction depended on whether the first 


random number was greater or less than pj. In like manner, the second random num¬ 
ber determined the direction along 9. 

Because step sizes were so small, the test-particle locations were not recorded 

every step. They were, instead, tallied at constant time increments At- Each time 

increment included Ati/(V, 0) steps. A count was kept of the number of times a test 

particle was found in each of the velocity zones (boxes). A large enough number of 

test particles A were walked to determine a velocity distribution. The time incre- 

max 

ments were selected so that usually no less than 10 nor more than 100 steps were taken 
between tallys. The number of Ar’s required to reach the loss cone was recorded to 
determine time of confinement and thus loss rate (appendix E). 


After each set of ste^ size and probability matrices were deter¬ 

mined with the test-particle distribution used as the field distribution. The elements of 
the tally matrix were reset to zero and the process repeated, if desired, for the i + 1 
step of the iteration. 

Numerous computations were made to explore the parameters of impact distance 

cutoff sample size A , and iteration index for suitable field convergence 

max 

I - I . External conditions such as mirror ratio and injection distribution were 
max 

selected either for reduction of computer time or for comparison of random-walk re¬ 
sults with cases having analytical solutions. 
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RESULTS AND DISCUSSION 


Influence of Impact Parameter 


The scattering angles, thus step sizes, descriptive of Coulomb encounters are, in 
general, so small that even to evaluate the effect of /3 with a reasonable amount of 
computii^ time, required selection of very short walks. Short walks were obtained by 
using low mirror ratios (B /B w 1). The initial condition for this study was set to 
simulate injection of monoenergetic particles normal to the magnetic field 0 = ti / 2 . 
Computations were carried out only for Z = 1 and a Maxwellian field distribution. 

A generalized parameter descriptive of the rate at which particles first reach pre¬ 
scribed conical boundaries in velocity space is plotted in figure 6. These boundaries 
(0 and 7T - 0 ) were selected so that no set of (either 100 or 1000) particle walks took 
more than 15 minutes of computing time. At least three values were used for each 
boundary condition. The mirror ratios (eq. (17)) became extremely close to one as 
bmax approached values representative of Debye lengths; in fact so low that these re¬ 
sults have little direct application to mirror systems of interest. These results may be 
viewed as determining the average length of walk required to first reach a certain conical 
boundary in a Maxwellian field. 

The regions of velocity space covered by the walks were so small that the step sizes 
and probabilities could be assumed constant over the length of the walks. This, in turn, 

Distance between start 
and finish of walk, 

Iwa -0^1, 

rad 

0.147 
.180xl0“J 
. 208x10"^ 

. 233x10"^ 

.294xl0"f 
.390x10"° 

—- Fokker-Planck solution 

{see appendix F) 

Open symbols denote 1000-walk sampling size 
Solid symbols denote 100-walk sampling size 


^ ° Q°o 0^8 ^ ^ Jl 

A 


2 3 4 5 6 7 8 

Logarithm of impact-parameter ratio, log (3 

Figure 6. - Effect of impact-para meter cut-off distance on rate of reach¬ 
ing prescribed 0 distance away from starting point at 7r/2. Singly 
charged particles of mass number 2; ratio of reference energy to 
field temperature, 15. 
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permitted an analytical solution of the Fokker-Planck equation (appendix F) which is 
shown by the line on figure 6. 

This analysis gave the result that the loss-rate parameter 

hT3/2(I 

log p 

was equal to a constant. This was verified over a range of /3 (from 10^'^ to 
by the computer results to within a spread of loss parameter of about 10 percent. 

Some typical test-particle distributions are shown in figure 7. The marginal dis¬ 
tribution in V 




vV„ 




f(0, v)sin B do - 



F{e, V)d9 


is a delta function in the analysis. It is sharply peaked in the computer results, nearly 
all of the points plotted lying inside 0. 998 :< v/v^ :< 1. 002. The marginal distribution in 
6 


Fg{e) = 


sin 9 


a 



f(9, v)v^ 


dv = 



F(0, V)dV 


for both the theory (appendix F) and the computer experiment can be represented by 
straight lines and agree well with each other. 

The theoretical distribution from appendix F is 





when 9 <0 < — 

^ 2 







when - < d < 7T - 9 
2 ^ 

J 


(18) 
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Marginal distribution in V, Fy^V) 



Angular distance from polar axis, 0 Velocity ratio, V 


(a) Confined distribution. 



Velocity ratio. V 


(b) Loss-cone distribution. 

Figure 7. - Test-particle distribution for short walks and constant coefficients. Ratio of reference energy to field 
temperature, 15; sampling size, 1000; impact-parameter ratio, 10^; distance between start and finish of walk. 
0.00208 radian. 
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The theory and computer experiment appear to be in satisfactory agreement for pur¬ 
poses herein and encouraged applyir^ the random-walk method to more complicated and 
practical models. 


Comparison With the First-Order Solution of Budker 

Analyt ical solution. - Among the very few analytical solutions of the Fokker-Planck 
equation applied to a magnetic mirror system is that of Budker (see refs. 5 and 19). 

This analysis makes the usual simplifyii^ assumptions of appendix B plus the assumption 
of separability of the test-particle distribution functions; that is, it is assumed that 

f(9,v) = f 0 ( 0 )yv) (19) 

inside the mirror system, and that the injection distribution (source) function is 


s( 0 ,v) = Sq{6)s^{v) 

-mv^/2kT. 

Solutions are then sought where fy(v) is of the form e . For injection of 

particles perpendicular to the B field the solution of reference 5 is 


-mv^/2kT^ 

f( 0 , v) = Ce 



-mv‘^/2kT. 9 

s{9, v) = Ce " <(A0) ) cos 9^ 5 


tl) 


( 20 ) 


( 21 ) 


where the constant C can be used for normalizing. Results were reduced to this form 
by approximation of an integration over 9 by the mean value theorem. The end-loss 
rate becomes 


^ 3 ln('^ + l) - Vi 

(47re )^(kT)^/2 Aj_\ 

« \sin 9j 


L 
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or 



In this solution the Fokker-Planck coefficients were based on a Maxwellian field 
distribution. It is thus like a first-order solution. 

No energy gain or loss from the system was considered other than that added by the 
injection process or carried away with the escape of particles through the mirros. Thus 
the mean energy of the particles injected should equal the mean energy of the par¬ 

ticles escaping for a steady-state solution. Int^ration of the product of equation (21) 
times 1/2 mv , and then division by ri gives the mean energy per test particle to be 

E. - = 0. 43 kTi. 

Thus the test particles are injected at a low energy to fill in the deficiency in velocity 
space due to the predominant escape of low-energy ions. 

Random-walk solu tion. - Some preliminary insight into the random walk can be 
gained from inspection of the step sizes and probabilities of taking steps in the various 
directions. Step-size parameters for the V and 6 components and for a Maxwellian 
field distribution are plotted in figure 5. Values ur^d in a typical set of walks are given 
in a matrix form in the computer listings presented in appendix E. 

For a given field distribution, path lines, as shown schematically on figure 1(b), 

can be envisioned. Use of the step sizes V<(^v)2 )/v and ) /v from figure 5 

permit study of the path lines for a Maxwellian field distribution. The magnitude of the 
slopes of such path lines increases with an increase of V because of the pronounced 

decrease of V(A0)^ with V. It is apparent that the path line pattern is most conducive 
to the escape of particles at low V. It is in turn relatively easy for a particle in the high 
slope region to go to still higher velocities where escape is more difficult 

The bias terms are equally important in the test-particle behavior. Enterii^ into 
the bias term (pj - qj) in the V direction (eq. (4)) is a dynamic friction term AV. 

This term (discussed in ref. 7) is negative and tends to retard the test particle to zero 
mean velocity. Opposed to influence of AV are the positive valued expressions 
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X (i9)2 (A,p)2 - i 

2 2 1/ 9V 

which tend to diffuse the particles to the thinly populated regions of space at high V. In 

2 

spherical coordinates the available volume in velocity space varies as V , making a 
large amount of room for particles at high V. 

Similar terms appear in the (p 2 - <12^ bias, which occurs in the 0 direction. All 

terms except the last (sin 9 cos 0/2)(A(p)^, are zero for a Maxwellian field distribution. 
For higher order solutions (I > 1) the AO and -A9 AV/V terms tend to move test par¬ 
ticles toward 9 = n/2. These terms thus help confine the particles. Opposed to this 

are the (AO)^)/dO) and [(sin 0 cos 0)/2] (A^)^ terms that make it easier for the 

particle to escape. 

The initial locations of the test particles were determined by selection of random 

numbers from a distribution to simulate the injection velocity pattern of equation (20) as 
shown in appendix G. The walks were always started at the 0 = 7r/2. 

Preliminary computer results verified that the test particles which were injected 
at low energies (according to eq. (21)) had a good chance of escaping with a relatively 
low number of steps. If, however, they remained inside the mirror system and by 
chance diffused to the high velocity region, their rate of escape was much reduced, even 
though the collision frequency was increased at high V. The test particles rush to such 
higher velocity zones at the expense of energy given up by the field particles. The field 
temperature was thus corrected for the energy exchar^ed with the appropriate test par¬ 
ticle after each walk (appendix E). 

Comparison of results. - Test-particle marginal distributions from the random- 
walk procedure are compared with the analytical results of reference 5 in figure 8. In 
general, the test particles distributed inside the mirror system in much the same pat¬ 
tern as predicted by reference 5 (figs. 8(a) and (b)) for the first-order results (i = 1). 

As the iteration process was continued, the test particles distributed in a slightly dif¬ 
ferent pattern. The peak of the distribution shifted to a slightly higher V. Beyond a I 
of 2, however, changes were within a 10- to 20-percent scatter of the data for the usual 
sample size particles. No attempt was made to optimize sample size 

against 13. The parameter /3 was usually selected for completion of times 

walks in 20 minutes of computer time. 

max 

Starting the Z = 1 calculation with the field-particle distribution in the 9 direction 
cutoff at 0 and ir - 0 appeared to give results as good as starting with a full Max- 
wellian distribution and completing the Z = 2 step. By so doing the final number of 
steps in the iteration process (^ could be reduced by one. 

The assumption of separability (eq. (19)) of the distribution function is often used to 
simplify analytical approaches (refs. 3, 5, and 19). This appears justified for the typi- 
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Generalized marginal distribution in 0,(1 - ?Q^l 7 r)^Q{Q) Marginal distribution in 
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in 0 direction 
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1 
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1.2 
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2 

Cutoff Maxwellian 
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3 
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- Maxwellian and ref. 5 

-l/2sin [(0-9(.)/(l-2 0j7i)] 
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□ 
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6 .8 1.0 1.2 1.4 1 

Velocity ratio, V 

(a) Test-particle distribution in V direction. 






(0 - 0(,}/(l - 

(b) Test-particle distribution in 0 direction. 

Figure 8, - Comparison of random-walk results with analysis of reference 5. Ratio of reference 
energy to field temperature, 1.5; sampling size, 500; atomic number, 1; impact-parameter 
ratio. 





(d) Distribution of particles as they first enter loss cone. 
Figure 8. - Concluded. 


cal set of data on figure 9. The open data points were obtained from the test-particle 
distributions of the random walk. These are a normalization of the KTJK(J, K) matrix 
(appendix E). The solid symbols denote the product of the random-walk marginal dis¬ 
tributions in 0 multiplied by the random-walk distributions in V. The lines are the 
analytical results of reference 5. The shift of the peaks of the computer results to values 
slightly higher than the analytical results is not due to the separability assumption. Good 
agreement is shown between the two studies for the outlet (loss cone) as well as for the 
inlet or injection distribution (figs. 8(c) and (d)). Thus the injection distributions of 
reference 5 appear quite satisfactory for a steady-state solution. 

A comparison of the random-walk loss rates with the results of reference 5 is 
shown in figure 10. Agreement is close at a mirror ratio of 1. 2. At higher mirror 
ratios, however, the loss-cone boundaries are extended enough that the chance of a par¬ 
ticle walking to high V fields before escape is considerably improved, accounting in a 
large part for the reduction of loss rate with increase of B /B^. 

The rather steep decrease of loss rate with mirror ratio obtained from the random 

walks quite closely follows the 1/1<:^(B /B ) trend predicted in references 3 and 13. 

The results of reference 5 differs by the inclusion of cos 9 in equation (22). 

c 


29 


[ 



F(e,V)or Fg(0)Fy(V) 


0 .2 .4 .6 .8 1.0 1.2 1.4 1.6 

8 - QJil - 2QJ7T) 

(a) Constant velocity ratio locations. 



Figure 9. - Study of separability of distribution. Sampling size, 500; ratio of reference energy to 
field temperature, 1.5; atomic number, 1; mirror ratio, 1.5; impact parameter ratio, 10^'^; 
iteration step number, 1. 




o 


Random walk with full Maxwellian in 
direction 



3.3 2.0 1.5 1.2 

Mirror ratio 

Figure 10. - Comparison of random-walk loss rates with the ap¬ 
proximate analytical solution of reference 5. 


The correction to the source term consisting of higher order terms in equa¬ 
tion (8) was evaluated numerically for a Maxwellian field distribution. For this case, 
equation (8) reduces to 


^corrected ^random' walk 


0. 003x10"^° log ^ 

mV2T3/2 


This correction, independent of mirror ratio is negligible in the random walk results on 
figure 11. 

Computer time expenditure. - Each complete set of A walks and I steps in 
the iterative procedure required less than 20 minutes of computer execution time. It 
appears that this random-walk procedure could be extended to include special or time 
coordinates, or both, and work more general problems with reasonable computer time 
expenditure. 


CONCLUSIONS 

From a study of computer simulation of random walks of charged particles through 
ensembles of field particles inside magnetic mirror systems, and by comparison with 
results obtained from the Fokker-Planck equation for cases that could be solved analyt¬ 
ically, the followii^ conclusions were reached: 

1. Computer time could be greatly reduced with no noticable reduction in accuracy 
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by a selective sampling of random numbers from the low impact-parameter range of 
the binary collisions used herein. 

2. Test-particle distributions were insensitive to the assumed field-particle dis¬ 
tributions. 

3. The assumption of separability of the velocity distribution appears suitable for 
application to mirror system analyses. 

4. End losses determined by the random-walk method were in close agreement with 
an approximate analytical solution of the Fokker-Planck equation at low mirror ratios. 
At higher mirror ratios the loss-cone boundaries increase permittii^ more test par¬ 
ticles to reach the higher regions of velocity space where escape is much more difficult. 
The loss rate varies inversely as the logarithm of the reciprocal of the mirror ratio. 

5. Each complete calculation, for a given set of initial and boundary conditions (in¬ 
jection velocity distribution and mirror ratio), was completed within 20 minutes of com¬ 
puter time. The study was limited to steady-state problems in velocity space. It ap¬ 
pears that this random-walk procedure could be extended to include spatial or time co¬ 
ordinates (or both) with reasonable computer time. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, October 20, 1969, 

120-27. 
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APPENDIX A 


SYMBOLS 


[The International System of Units (SI) will be used throughout with the exception of 
temperature T, reported in keV, and the corresponding Boltzmann constant k, re¬ 
ported in J/keV. ] 


A 

a 


B 

b 

C 

E 



e 


F(0,V) 


F,(a) 


Fy^V) 


f 

g 

H 

h 

k 

I 

m 


mass number 

determinant of elements in metric tensor 
elements of metric tensor 


elements of cofactor of metric tensor 

magnetic field 

impact parameter 

normalization constant 

energy 

mv^/2k 

electric unit of charge 

normalized probability density times scale factors, 


sin e 


n, 


f(v,0) 


normalized marginal distribution in 



F(0, V)dv 


normalized marginal distribution in 



F(0, V)d0 


probability density 


d^ 
b 

Heaviside unit function 
Roseribluth potential = 

Boltzmann constant = 1. 6x10’^® j/keV 
step number in iteration process 
mass 


Rosenbluth potential = 
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N 

n 

* 

n 


Pi 



gliv 


9 

s 


r, R 
T 

tM 

t 

u 

u(|j, ^2’ 

V 

r 

V 

V 
Z 


r 

6 

e 


number of dimensions 
number density 

number of particles lost per unit time per unit volume of coordi¬ 
nate space 

th 

probability of 1 step in the positive direction along the i coordi¬ 
nate 

th 

probability of 1 step in the negative direction along the i coor¬ 
dinate 

general curvilinear coordinate 

r"^<Aq^ Aq^> 

number of particles injected per imit time per unit voliune of phase 
space 

random numbers 
temperature 
r"^<Aq^) 
time 

magnitude of relative velocity between a field and test particle 
probability of particle being located at grid location 

^ 1 ’ ^ 2 ’ • ■ 

v/Vq 

random number 
velocity 

some arbitrary collisionally dependent quality 
atomic number 
impact-parameter ratio, 

(Zg)^ In 

interaction parameter, =- 

^22 
4776 ^m 
o 

delta function, also used to denote finite difference 
azimuthal angle between orbital and fundamental collision planes 
capacitivity of vacuum 
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7] number of components changed per encounter 
6 angular distance from polar axis 

angle between test and field-particle velocity vectors 
u collision frequency 

^ arbitrary coordinate in phase space 

a differential scattering cross section 

cp azimuthal coordinate 

Q solid scattering angle 

Subscripts: 

b field particles 

c collisional value, also used to denote loss cone or value at mirror location 

e electron 

inj inj ection 

o reference value 

Special symbols: 

average per collision or per step 
() time average 



APPENDIX B 


TRANSFORMATION OF THE FOKKER-PLANCK EQUATION FROM 
RECTANGULAR CARTESIAN TO SPHERICAL COORDINATES 

The Fokker-Planck equation can be written in rectangular Cartesian coordinates as 

2 

= - _L (f<Av^>) +i -^- (f<Av^ Av^>) 

'^^ 4 ; ^ 9v^ dv^ 

Superscripts are used in this appendix to denote component direction, leaving subscripts 
to denote covariant differentiation. Summation convention of repeated indices will be 
used. The subscript c on 0f/9t means changes of f due to collisions. 

Following the procedure of appendix IV of reference 11, this equation can be written 
in covariant form as 


Vaty ’ 2 


IJ.V 


12 3 

which is valid for any set of curvilinear coordinates, q , q , q . 
The quantity F is an interaction parameter. 


tM = ^ = <Aq^>r 


-1 


(Bl) 


and 


= (Aq^ Aq^> r"^ (B2) 

The terms are dispersion coefficients, and TT^ is related to dynamic fric¬ 

tion. 

12 3 

For spherical coordinates q = v, q =9, and q = (p. The elements of the 
metric tensor are 
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agg = sin^0 


a = 0 

IJ-v 


if IX ^ V 


and the elements of its conjugate are 


all=l 


22 1 
a = — 


„33 1 

a =- 


2 . 2n 

V sm e 


a^'^ = 0 


if IX * V 


By the procedure of covariant differentiation 


(fT^) 


M ^" {a*;.}^ 


r“ 

where the Christoffel symbol { ^ I with repeated indices equals ——— and a is 

va a„ O' 

3q 

the determinant of the matrix of elements a . For spherical coordinates 


4 9 

a = sin^0 


Thus, 


r(fT^) = Ji--L (v^f(Av>) +— - -^(f sin 0(A0)) +-^ (f(A<p>) 

’ ^ 2 0v ' sin d dd dcp 


In like manner 


r(fs^^), 


1 

va V ^ / 


3q/^ aq*^ 
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and for spherical coordinates 



WiK") • 



sin e ^ 

2 \ 


0a 9a ^ \ 
0q^ 0q^ / 


— fv^ sin 9 f(vS^^ + V sin^0S^^)l + A [sin 6 f(2vS^^ - sin 6 cos 0 
dv^ -' 00 ^ 


+ r2fv(sin 0 + V cos 0 

0(/? L 


where use is made of the symmetry of S^. Finally the Fokker-Planck equation in 
spherical coordinates becomes 



— — (v^f <Av)) 

v2 av 


-^A(f sin9<A0)) --L(f<Ay))^J_g^ f((Av)^) ^ 1 a^f sin 9((A9)2) 

sin 9 39 3<p ^^2 g^2 2 sin 9 gg2 


+ - \ —I - 1 -(v^ sin 0f(Av A0)) - (v^f(Av A<;^?>) +—^ -? 

9v 00 ^2 0v d(p sin 0 00 d<p 


dcp' 


v^sin 0 


(sin 0f(A0 A(p)) 


— — pf(((A9)2)+sin20((A(p)2))] +-^- 

2^2 0v ^ 2v sin 0 


A[f sin 0 ^2(A0 Av) - v sin 0 cos 0((A^)^)^j 


- f Av A(p) 

V 0<^ L \ 


^ ^ <A9 Acp}) 

sin 0 /- 
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APPENDIX C 


PHYSICAL MODEL OF MAGNETIC MIRROR SYSTEM 

The physical model to which this analysis is applied is the usual simplified one, 
well described in many references; such as, for example, pages 186 to 189 of refer¬ 
ence 5, pages 3 to 5 of reference 13, pages 2 and 3 of reference 3, and pages 7 to 9 of 
reference 17. In essence, the magnetic mirror system provides a long cylindrical re¬ 
gion of uniform magnetic field (fig. 3). The mirror regions are assumed to be rela¬ 
tively short so that essentially all the Coulomb scattering occurs in the central re¬ 
gion. Only binary collisions and classical particles will be considered. 

Once a particle enters the loss cone in velocity space, no matter where it is in 
physical or coordinate space, it is assumed lost from the system. Only scatterii^ of 
ions off other ions is considered. Thus electrons serve only as a neutralizing and 
shielding bacl^round of particles. For study herein, no macroscopic electric field will 
be accounted for except possibly as a boundary condition in the manner of reference 20. 

The radius of the confining field must obviously be much larger than the ion cyclo¬ 
tron radius. This can be accomplished by suitable magnitude of the B field. The 
azimuthal symmetry of the B field makes v x B • V^f = 0 as shown in appendix A-1 
of reference 13. 

These assumptions culminate in the elimination of the gradient terms, in velocity 
as well as in coordinate space, from the Boltzmann equation. The Boltzmann equation 
thus reduces to 


at 



+ s 


where (0f/at) is the change of f due to collisions, usually determined by the Fokker- 

V./ 

Planck equation, and s represents a sink density in velocity space, uniform in coordi¬ 
nate space. 
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APPENDIX D 


FOKKER-PLANCK AND RANDOM-WALK COEFFICIENTS 
IN SPHERICAL COORDINATES 


The Fokker-Planck coefficients are time averages of various first- and second- 
order increments of velocity due to collisions. The random-walk process, in turn, re¬ 
quires the average per collision of these same velocity increments. The two types of 
averages, denoted by the angular brackets and bar notations, respectively, differ by 
the collision frequency v. Thus only one set of coefficients and the collision frequency 
need be derived for the physical conditions of interest. 

The time rate of change of some collisionally dependent quantity Y averaged over 
all scattering angles and all velocities of the scattering (field) particles is 


<Y> =A(vi,)/Yua( u, 


(Dl) 


where the subscript b refers to field particles, and a is the differential cross section 
of scattering over the solid angle f2. The magnitude of the relative velocity between a 
test-particle of velocity v and a field particle of velocity ^ is 





cos 

r 


(D2) 


where is the angle between v and v^^. The distribution function of field particles 
is denoted by f^ and is normalized so that 

Using Rosenbluth potential (ref, 11) defined as 


= /f(v,)idv. 


and 


g = 
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■fVi 

and letting Y equal the component of Av yields (ref. 11) 


(Av^> = r 


8h 

9v^ 


i = 1, 2, 3, 


the covariant form of which is given by equation (Bl) 


<Aq^> = ra^‘^(h 


For Coulomb encounters 


4 4 

„ _ Z^e^ In i3 

.22 

47re m 

where /3 is the ratio of the maximum impact parameter b„„„ to minimum, b„. . 

• • XIxaX IXxlTl 

When Y is set equal to the product Av^ Av^, it follows that (ref. 11) 

.2 

(Av^ Av^> = r-S— i, j = 1, 2, 3 

• « 

0v^ 9v^ 

and the covariant form is given by equation (B2) as 

<Aq^ Aq*^) = ra^‘^a'^'^(g^ 

When Y = 1, it is apparent from equation (Dl) that the collision frequency is ob¬ 
tained. Using the well known scattering relation involvii^ solid angle fi, impact pa¬ 
rameter b, and azimuthal angle e 


a’(n)dn = b db de 


in (Dl) gives 


V 


irfb 


2 .2 

- b . K 
max min/® 


) 


which is an invariant and independent of the coordinate system. In spherical coordi¬ 
nates, the coefficients of main interest are 



<Av> = r — 
9 v 


<A0) =L.^ 
..2 dd 


(A<p) - 


r 3h 

2 . 2^ d<p 

V sin 0 ^ 


(Ad Av)=i:fj!g_. iii' 

2\00 av V 00 , 


<(Av)^) = r-^-i 

0v2 


<(a0)2> + 

vAa02 


((A(p)^} =--/A_i_ + V sin20 ^ + sin 0 cos 0 ^ 


4 - 4 . 1-2 

V sin d\d(p 


0 V 


00i 


<AvAy>._r 

v2 sin20 \ ^ 


(Ad Acp) =_jg^ 

sin20 \ d dcp) 

These expressions can be reduced one step further without specifying f^ since differ¬ 
entiation is with respect to the coordinates of the test particles. 

Azimuthal symm etry. - For aximuthal symmetry d/dcp is zero and thus (A^), 

(Av A^) , (A0 la^cp) are zero. The remaining expressions are reduced as follows: Using 

cos = cos 6 cos + sin 6 sin 0^ cos{cp - cp^^) 
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and differentiating under the integral sign gives 


<Av> = -2r 



V - cos dj. 


- 2vv^j cos 


3/2 


f(Vb)dv^, 


(9a) 


<A0> = - 


2 p y w^jj^sin 9 cos 0^ - cos 9 sin 0^^ cos{(p - (p^)j ^ — 


(v^ + - 2w^j cos 


3/2 


(9b) 


(A0 Av) = — 
2 


p cos B^) j^sin 0 cos 0^ - cos 0 sin 0^^ cos(<p - 


(v2 + v2 - 2vv^ cos B^) 


3/2 


f(%)dVb 

(9c) 


<(Av)^> = r 



(V - cos ap 


- 2vv^ cos ^ ■ ^^b ^r) 


3/2 


f(vpdv^ 


(9d) 


((Ae)2> =L. 



r 


2 2 r 

V sin 9 cos 0^ - cos 9 sin 6-^ cos 


w - <Pb)] 


(v^ + - 2w^ cos 


' 2 ..2 


V + - 2wj^ cos 9j, 


3/2 


f(Vb)dVjj (9e) 


<(A<p) 2) = r 


sin^0 



V - vv^ sin 0^ cos(<p - <py^ 


sin 0 


V + - 2W|j cos B^ 


1/2 


f(Vb)dv^ 


(9f) 


4r)3^ / m 
V = —c— 


9 In ;8 V2kT^ 


■) / P + 'b - cos S. 


«Vb)dVb 


(9g) 


where b^^^- was set equal to the classical distance of closest approach 
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(Ze) 


_ (Ze) 


mm 


o(l"») 


47re 


6;r€okTb 


The average value of these quantities per collision, in dimensionless form, are 



A 0 AV = 


VV, (V - Vt^ cos 5 ) sin 6 cos d. - cos S sin 9-. cos((p - ^ 

— 2 -E- T L —.-P_. . ... oj f{vJdVK 

3/2 ^ 


<A 0 AV) 

9 In /3 j 

^Tbf J 

(v 2 + V^ - 2 VVjj cos 3 J 

V 

%v’ 


/Vv 2 + v 2 - 2 VVj, cos 3 , 


(D4c) 


(AV)2 = = 9 togPb 



4 /32 


(V - Vj, cos 3 j.)' 


^V 2 + v 2 -, 2 VVj, cos 3 ^ (v 2 ^ - 2 VVb cos 


f(Vb)dVb 


/V^ 


Vb - 2 W^ cos 3 ^ 


(D4d) 


(A0)‘ = 


2 _ <(A9)^> 


_ 9 In ^ / b \ J 


V^V^pin 0 sin 0 jj - cos 0 cos 0 jj cos(<i!; -‘?’b)] 


/v 2 + v^ - 2 VVjj cos 


(v2. 


Vb - 2^b \) 


3/2 


f(Vb)dVb 


4 v V 




V‘ - 2 VVb cos Sj. f{V,,)dVi, 


(D4e) 
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In a 


(A(p)^ = 1 = j 


,EL 



V - VV, 


Q 

/sin 0 sin 0, + cos 0 sin 0^ 

-- -- -- --^ cos(<p - (P^,) 

’V sin 9 /_^dV 




^ v‘^/3^ sin^e 


-f - 2VV^ cos 3^ 
y-H - 2VV^, cos flj. fCVdV^ 


f(Vb)dVb 


(D4f) 


and 


V = 


4rj3 
9 In j3 


2 / \^/^/e' V/2 


f—) 

\2kT / 


/V 


- 2VV^j cos 


(D4g) 


where 


V 


Vb = 


E = i mv^ = kE' j3 » 1, 


o o 


(D5) 


Maxwellian distribution . - A popular simplifying assumption is that the field par¬ 
ticles are distributed in a Maxwellian spherical distribution. This a case of the ’’diffu¬ 
sion approximation” on page 175 of reference 5 and corresponds to the use of ’’effective 
Rosenbluth potentials” as on page 15 of reference 19. With this assumption 




= n. 


m 


3/2 


-mvn 


72kT. 


27rkT. 


sin 9 . dft 


dcPb^Vb 


( 10 ) 


For a spherical distribution, 5 can be replaced by 0^. Integrating between the limits 

0 < < «> 

0 < 0^ < ff 


and 


0 :< < 27r 
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gives 


g = g(v) = njj 


, If 7rm 


2kT^j -mv^/2kT^j 

e 


kT 

+ IV + —lerf f ■* /-51L v 


mv/ \ If 2kT 


(Vi 


The remaining potential h can be obtained from 


-..2 

^ V \I2kTu 


Since both g and h are independent of d and <p, it is apparent from (Bl) and 
(B2) that 

(Ad) = (A<p) = (Ad Av> s 0 

The remaining terms of interest for this field distribution are 


(Av) = 


2n^r 


V 


i 


m 


27rkTi 


ve 


-mv'^/2kT^^ 


“ erf 


(V^'' 

V\2kTj, ! 


(11a) 


kT 

<(Av) 2> = _ _b <av> 


mv 


(lib) 


<(A0)2) 


1 

f Tim V 


^ -mv^/2kTjj / 


kT, 


mv 


erf 


s 


m 


2kT, 


(11c) 


<(A<^)2> =_1_<(a0)2> 


sin^0 


(lid) 


and 


u = 


9 In (3 \2kT, 


vr 


l2kT^ 1 -mv‘^/2kTv, / kT 

—^ie ^+(1 + 


Tim V 


mv 


.2/ V\2kT^ 


(lie) 
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2 

These expressions for (Av> and {(Av) ) agree with those of reference 5 and with 
those of Chandrasekhar (pp. 73 to 75 of ref. 12). These coefficients are called paral¬ 
lel components since they are changes in v in a direction parallel to v. Agreement 
with these references is also obtained for the perpendicular component of diffusion 
which is 


<(AVj^)^> =v^|^<(A0)^> + sin^0<(A^)^>] 
Remainii^ quantities of interest are 




+ 9 W-” 

2 / VEo. 3 

yO 




e o ^ - erflV 




- 

Ml 

f^b 

v \ 




2V 




n 


2 


(AV)^ = - 


2VE’ 


AV 


(15a) 


(15b) 


(A0)^ = 


<(A0)2) 


9 In^Ab 
4 ^4^2l,E'. 




o 


. .1 - 


T 
2VE 


m'J \ Vt^/ 


-V^E^T^ 
e + II + 


b 


2VE’ 


o> 




(A<p)2 = _J_ (A0)2 
sin^0 


(15c) 


(15d) 


and 


u = 


0. 518x10"^^ n^zVv 

T3/2ml/2 




1 + 


2V^e! 


srf V' 


O' 




(15e) 


where B » 1. 
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APPENDIX E 


DESCRIPTION OF COMPUTING PROCEDURE 

A flow chart representing a typical calculation procedure is given in figure 11. The 
corresponding FORTRAN program and some typical results are included at the end of this 
appendix. This is preceded by a key to the computer words. A brief description of the 
prc^ram to perform the random walk studies is given in the following discussion. 

The distribution of field particles is represented by KTJK(J, K) matrices. The 
indices J and K represent partitions in V and d, respectively. In most of the com¬ 
putations this was an 11 by 13 matrix, thus dividing the velocity space of most impor¬ 
tance into 143 zones. The initial (LL=1) field distribution was read from data cards. 

The velocity field to which the test particles had access was also divided into zones. 
These zones were identified by subscripts M and L. Step sizes, DVSQAV(M, L) and 
DXSQAV(M, L), and probabilities, TWOP(M, L) and TWOR(M, L), were computed over 
the ranges of M and L for the field distribution specified by the KTJK matrix. Inte¬ 
gration of equations (9) with respect to <p was by Gaussian quadratures. Integrations 
with respect to V and 9, weighted by v sin 9 f(v, 9) as in equations (9), were replaced 
by summations over the J and K indices with weighting by KTJK(J, K). 

Because step sizes are so small, the test-particle locations were only recorded 
every MMAX(M, L) steps. To obtain steady-state distributions, it is desired to locate 
the test particle after each small constant increment of time. At. The number of steps 
per constant increment of time is 


RMAX(M, L) = (RMAX1)(j/(M, L)) 

v{Y2=\, X2=ff/2) 

where RMAX is the real number counterpart of integer MMAX. The value of RMAXl 
is specified at the beginning of the program and should be judiciously selected to keep 
computing time down and yet obtain accuracy. 

An alternate procedure to the use of velocity zones and matrix descriptions when 
LL=1 was to specify a Maxwellian distribution of field particles and use equations (15). 
In this way the step sizes, probabilities, and number of steps per time increment are 
continuous functions of test-particle location and are integrated over a continuous field 
distribution. A comparison of this procedure with the matrix procedure served as a 
check on the required number of velocity zones for sufficient accuracy in the matrix 
representation. Excellent agreement was obtained by use of the 11 by 13 matrices. 
Suitable values of RMAXl varied from 10 to 10 depending primarily on 

The test particles were labeled by A (used as an integer). It was found that the 
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walks of about 500 test particles (AMAX=500) was required for a representative sample. 

The initial conditions for the random walks depended on the injection distribution 
being simulated. For the example on the flow chart, V2 was always set at 1. 0 and X2 
at 7r/2, representative of monoenergetic injection and normal to B field. This initial 
condition was used in the short walk calculations. For simulating the injection distri¬ 
bution of reference 5, the procedure of appendix G was used. 

The random-walk proper comprised a small part of the programming effort. In 
the two-dimensional space of interest it required the random selection of 2 numbers, 
RNR and RNS. The call, RAND (RNl), for the first random number, RNl, must be 
preceded by SAND (RNl). This sets up addresses in the congruence type of random- 
number generator subroutine of the Lewis computer library. This is a pseudo-random 
number generator inasmuch as each time SAND (RNl) is called the same sequence of 
random numbers is started again. Comparison of RNR and RNS with TWOP and TWOR 
determines in which directions the steps are taken. Step sizes and probabilities are 
held constant over MMAX(M, L) steps. After MMAX steps are taken the new location of 
the test particle is determined and tallied. Also after MMAX steps the energy transfer 
between the field and the test particle is determined and summed to that of the preced¬ 
ing N groups of MMAX steps. 

If the test particle has not yet reached the loss cone, it is tallied according to the 
procedure on the lower right of the flow chart. It is simply located in one of the veloc¬ 
ity space zones and credited to the corresponding element of a KTJK matrix. A cer¬ 
tain probability exists for a particle to wander about indefinitely inside the mirror sys¬ 
tem without being lost. A walk length of 5000 MMAX steps was therefore set as a limit, 
after which the walk would be terminated and the next particle selected. The number 
of such walks was labeled IBAD. 

When a test particle reaches a loss-cone boundary, the flow in the chart moves to 
the left. If it is the AMAX’^^ particle, the tallying procedure in the lower left is fol¬ 
lowed. The marginal distributions in V and d are then determined along with a 
counting of the total number of points (KTOT) tallied in the KTJK matrix. The sum of 
the values of N when the particles first reached a loss-cone boundary is printed out 
as SUM of NA. 

Denote the collision frequency per test particle at X2 = tt/2 and V2 = 1 by 
t>(7r/2, 1); and define QUCOR by 


■'(? 

QUCOR = —L - 

^b 



_m_ 

Mass of a deuteron 


Then by use of equation (9g) 
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CALL SAND RNl) 


READ EOPT, BMAX, RMAXl, AMAX, LIVIAX 
XMIN, XMAX, VHI, XHI, DELVHI, DELXHI 


READ KTJK(J,K) J» 1. 2.11: K = 1. 2 


LL« 1 


CALCULATE AND WRITE DVSQAV{M, L) 

DXSQAV(IV\.L), TW0P(M,L)TW0R(M.L), MMAX(M,L) 
ANDQUCOR. M» 1.11; L» 1, 2.13 


KMARGV(J) = KMARGX(K) = KTJK{J, K) = 0 
KTOT = SUMNA « IBAD » DVNSQ = 0 
TR = 1.0 


\ VINJ » V2 » 1.0, X2 » 7r/2, N = 1, IX = lY - 0 / 


A-A + 1 

























XHT-XHI 

SUMNA • SUMNA + FLOAT(M) 

No 


VOUT -V2 


VHT • VHI 

Tp _ 1.0 “ DvNSQ + FLuAi(a) 

“■ ■ * . V^MIN V X2^ XMAa^ ■ • 


XOUT-X2 


J-1 

SUMNA - aOATfAMAX) 





K-1 


KMARGVtJ}- KMARGV(J) + KTJK(J,K) 



Tally for marginal 
distribution in V 


I Tally for marginal 
distribution in 6 


Tally for field 
distribution 


Sum of points tallied 


WRITE KTJK(J.K), KIVIARGX(K) 

KMARGV(J), J = 1, 2 ,. . 11; K= 1, 2.13 

KTOT, QUCOR, SUMNA. AM 




K-K + 1 

XHT - XHT + OaXHI 


j-J + l 

VHT ' VHT + OaVHI 


iVlM= 1 
N-N + 11 


IBAD-IBAD + 1 
A -^A ■ 1 


I 


LL-LL + l 


Figure 11. - Flotv chart. 
























QUCOR = 0. 896x10" - 

KTJK(J,K) y „ 2V SK ■ cos X. 

KTOT 1 » D i 

J, K 

i 


for z = 1. The average number of Ar time increments per walk is SUMNA divided by 
AMAX. The average number of walks per unit time multiplied by 

^T^/2y4^^m/Mass of a deuteron is then equal to (QUCOR ' AMAX)/(SUMNA • RMAXl), 
and the loss rate parameter can be expressed as 

_ 0. 578xl0~^^ QUCOR • AMAX 
^2 SUMNA • RMAXl 

After each walk, the field temperature is corrected to account for the energy trans¬ 
ferred with the test particle. Assuming AMAX particles in the field ensemble results in 



_A_ 

SUMNA • AMAX 


SUMNA 

i: 


(V2'^ - 


VINJ^) 


N=1 


(E2) 


The correction is applied through the dependence of step sizes and probabilities on 
v'^m/2kTj^ as in equation (10). 

After A reaches AMAX the resulting test-particle distributions are used to calcu¬ 
late new step size and probability matrices, and the iteration process in LL is continued 
until LL reaches LMAX. 


FORTRAN SYMBOLS 


A 

test-particle number 


AM 

AMAX - IBAD 


AMAX 

size of test-particle sample 


BMAX 

impact parameter ratio 


BRATIO 

mirror ratio 


CK 

cos(X) 


CRELi 

cos 6 cos 0^ + sin 9 sin cos(7Txp 
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cx 


cos(X2) 

CX*CK 


CXK 

DELTAV(J, K) 

DELTAX(J, K) 

DELTHA 

DELVHH 

DELVHI 

DELVSQ 

DELXHI 


(AV)(J, K) 

A?(J, K) 

(A0)^ 

DELVHI* TR 

increment of V between successive J indices 
(AV)^ 

increment of Q between successive K indices 


DENOMl 53 KTJK(J, K)G]S1D{J, K) 

J, K 


SUMNA 


DVNSQ 

(V2^ - VINJ^) 


N=1 

DVSQAV 

1 —————— 

V(av)2 

DXSQAV 


EOPT 


FDV 

1.21/(SIB*EXP(Z*Z)) 

FDVi 

Ti/Ui when i = 1, 2, 3, 4 

FDVSi 

TT2Vui 

FDXi 

V*SRELi/Ui 

FDXSi 

1.0/FNUi - Y*SRELi**2/U2 

FFVi 

FNUi* FDVi* FDVSi 

FMIXi 

V*SRELi*FDVi 

FNDi 

SQRT(1.0+Y-2.0*V*SK*COS(7rXj)) 

FNUi 

SQRT(W+Y-2.0*V2*V*(CRELi)) 

FPSIi 

(Ti+V*CX*SRELi/Sx) /FNUi 

FRICH(M, L) 

2.25*S*NUM9/(W*NUM1) 

FTERM(M, L) 

-6. 75*S*NUM6/NUM1 
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FTHEi 
F (THETA) 
FV 
F(V) 

GDV 

GDVS 


V*SRELi*(V2-2.0*CRELi-3.0*V2*Y(SRELi/FNUi)**2) /Ui 
KMARGX(K) print out 
ERF(Z) -R 

KMARGV(J) print out 

Vw.FDVi(x.) 

i 

y^w.FDVSi(x-) 

i 


GDX 


2]w.FDXi(x.) 


GDXS 


GFV 


GMIX 


Vw.FDXSi(x.) 

i 

2]wjFFVi(x.) 

i 

23 w.FMIXi(xj) 

i 


GND 


J]w.FNDi(xp 

i 


GNU 


jy.Fmiix.) 

i 


GPSI 


^w.FPSIi(xj) 

i 


GTHE 

IBAD 

IX 

lY 


J]w.FTHEi(Xj) 

i 

number of walks discarded due to N reaching value 5000 
number of steps in positive 9 direction 
number of steps in positive V direction 
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KJ(J) 


KMARGV(J) 

KMARGX(K) 
KTJK(J, K) 

KTOT 

LJ(J) 

LL 

MMAX 

N 

NA 

NUMl, NUM2, 

QUCOR 

R, RNl, RNR, 

RAND 

RMAX 

RMAXl 

S 

SAND 

SIB 


4* V\ 

niimber of points in interval of injection distribution 
in V 

^KTJK(J, K) used for marginal distribution in V 
K 

^ KTJK(J, K) used for marginal distribution in 0 
J 

f f Vi 

tally of points in increment of V and increment 

of 6 for field distribution determination 

£ KTJK(J, K) 

i,k 

VVi 

number of points in increment of V in loss-cone 
distribution 

step number of iteration process 
number of steps between tally points 
number of groups of MMAX steps 

number of groups of MMAX steps when particle reaches 
loss-cone boundary 

. . ., NUM9 used in weighting functions by KTJK and summing over 

J and K 

collision frequency at V2 = 1. 0, X2 = 7r/2 
RNS random numbers selected from uniform distribution 

call code for random numbers 
real MMAX 

MMAX when V2 = 1.0 and X2 = 7r/2 
In 13 

initial call code to set up addresses in random number 
generator 

Vv®o 
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SK 

sin X 

SRELi 

sin d cos 6-^ - cos 9 sin 0^ cos{rrx^ 


AMAX 

SUMNA 

IT na 


A=1 

SX 

sin(X2) 

SXK 

sin(X2)sin(X) 

Ti 

V2-V*CRELi 

TR 

1 - AT^/Tjj, see eq. (E2) 

TTi 

V*SQRT(1.0-CRELi**2) 

TWOP 

2pj 

TWOR 

2p2 

Ui 

FNUi**3 

V 

magnitude of field-particle velocity 

V2 

magnitude of test-particle velocity 

VHH 

VHI*TR 

VHI 

initial location on V scale for tallying particles 

VHT 

initial location on V scale for tallying particles 

VINJ 

injection velocity 

VO 

Vl. 5 T|j/E^ arcsin(R^) 

VOUT 

V2 value of point being tallied 

VSQFV(J) 

KMARG(J) 

VSQFVT 

K]V1ARGV(J) 


J 

VSQKJ(J) 

V^KJ(J) 

VSQLJ(J) 

V^LJ(J) 

VSQKJT 

KJ{J) 


J 

VSQLJT 

Ev^ LJ(J) 


J 
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w 


(V2) 


2 


X 

X2 

XHI 

XHT 

XMIN and XMAX 
XOUT 


Y 

Z 


Subscripts: 

J 

K 

L 

M 

i 


weight factors in Gaussian quadrature integration procedure 
abscissa location in Gaussian quadrature integration procedure 

®b 

9 

initial location on 9 scale for tallying particles 
initial location on 0 scale for tallyii^ particles 
loss cone boundaries on 0 
X2 values of a point being tallied 


1.073 

index 

index 

index 

index 

index 


voVeJTt^ 

on V or V2 
on X or X2 
on V2 
on X2 

on Gaussian quadrature terms 
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FORTRAN PROGRAM 


C I^^ -NS ION KMAj<GV( ll)tKf^ARGK(13),KTJK(ll,13),KJ(ll),LJlll) 

Cjy _NS f ON PvSGAy ( 13 Lj CKSQ^VH LtXii tlhllP i l U U 1_» JHOAi 11, 13 )» 
L ll»13)v VSOFV(ll) 

CIMcNSIDN OtL[AV( LL,^iJ.),L:ELTAXlIl»13J,L)PSI ILL, ii) TCc<M(11,13 1 
IFRICH(11,13) ,CMIX(11, 1 3 ) 

INTAQFR A,AMAX 
CALL SAND(RNI) 

^^tAD I\ OF INITIAL field C I S TR I OUT I O'i ... .. 

REA,:(5,502){(KTJK(J,K),J=1,11),K=1,13) 

Ar lSi 5Q3J LKKAKGVCJ ) ,J= li 11) . 

RhAOIStSOA) (KhARGX(K) ,I<=1, 13) 

PEAL L5,50 3 )KT.UT 
502 FCR^AK 1114)' 

F.QRf-AT( 1115] „ .... 

:j04 FCR'^AK I 315) 

305 FCRi^ATIlil .... 

UR I rF( 6,625 ) 

. . >RITE(6,_62.‘^} IKi^tKTJKU.i; ).. J= L, 11) » Ki^ARGX (K ) . K= 1,1 3 ) 

URI(F(6,629)(KMARGV(J),J=1,11) 

_ _.kRI ... .. 

10 PF Al. (5,500 ) tCPT, LMAX, RMAXI , AjVAX, LMAK 
.If tLCP_T.EG.*OMlGu. TO 17Q... . . 

PEAl( 5,501)XMIM,XMAX,Vhl ,XhI,UELVH[,OELXHI 
. URI rE(.6^6CC)EUPT4-bMAX»RiX.AXl, AMAX,XMIN,XPAX . 
PRATIG=1.0/(Si:'MXiMIN) 

J?RI1JZ16*63U Viir.i..XHlji].aVJ:J.jiD£LXi±I.,JiR.AJJLi^ . . ^ _ 

CC30LL^1,L^AX 
TR=l1^P . 

S = ALnG( LMAX)/(FUPT*BMAX) «*2 

^ I TERATIj'NJ)Q_ LOUP FOR SELF C CN S I ST FNT„ F f EL D .01 SIR I BUT I OH 
CC29L=1 , 13 

C CAIC L LATIGN. ..Ot hJkP AfsH^PfOr AEat IJY NATRIPE S.^.. _ 

f:C?5N^=l, 1 1 

V2=>/Hl + (FLf^AT(M)-l.5)*DELVHI ^ . 

X2=XHI+(FLCAT(L)-1.5)»CELXHI 
U = VZ**'? 

CX = CGS(X2 ) 

SX = SIN(X2)_... ... .. ..... _ 

REAL NU^l , NUM2 f NUN3, MUN4 , NUM 5 , NJUM6 , NUM? , NUM 3 , WUM9 
. ^UY1=0.0 ... 

M9^2=0.0 

NUP3=0.0 . ... 

^UM/^=0.0 

_... . , . . _ ^ ... 

5=0.0 

^U^l/=0.0.. .... ... 

NU^^>^=0.0 

i=a^Q . . .... __ .... .. 

CF.\CM1 = 0.0 

_JNTEGRA„UOil. QVfil.PCLAj<.ANGtE__41vU. YELQCITY... B_Y QUADRATURES 
CC31K=1,13 







rc3?j=ifii . _ __ 

IF(LL-1)13,13,14 

13 V= C 0-200+ tFLCUTl J)-l^ 5)*Ce200 J *SQRT(1* 5/EOPT } 

>=0.241661+(FLOAT(K)-1.5)*0.241661 

GC TO 15 . . __ 

14 V=VHI+(FLOAT(J)-1.5)*CELVHI 

_ )( = Xh T-K FI PAT (KJ-1 »rFI XH T _ 

15 V=V**2 
rK=cnstx) 

SK=SIN(X) 

CXK=CX*CK 

SXK=SX*SK 

r. TMTFFRATinN DVFrt A7fMIJTHAL A^GLE BY GAUSSIAN QIJAnRATURES _ 

CRELl=CXK+SXK*CaS10.16343*3.1416) 

SREL1=SX*CK-CX*SK*C0S C 0.18343*3* 1416) 

FNU1=S0RT (W + Y--2.0*V2*V*CReLl ) 

PNOl = SQRTf ie-0+Y-2^0*V*SK*CCS (0^18343*3.1416) ) 

L1=FNU1**3 

T I=V?-V*rRFL 1 ________ 

TT1=V*S0RT(1.0-CREL1**2) 

FOXl=y*SRELl/Ul ^ 

FDV1=T1/U1 

J^DVSl=TTl**2yLil - _ 

FDXS1=1.0/FNU1-Y*SREL 1«*2/U1 

_ FF V l=F^!Ul *FDV1«FPVS1 __^_ 

FPSI1=(T1+V*CX*SREL1/SX)/FNUl 

FMTXl=V*SRFLl*FOVi * - - . 

FTHEl=V*SRELl*(V2-2.0*CRELl-3.0*V2*Y*(SRELl/FiNjUl )**2)/Ul 
CRFt2=C XK + SXK*C£}S (0^52553*3. 1416) 
SREL2=SX*CK-CX*SK*COS(0.52553*3.1416) 

_ FNU2=SQRT ( Q*V2* V*CREL2 ) _. . 

FND2=SQRT(1.O+Y-2.0*V*SK*COS(0.52553*3.1416)) 

L25sFN1I2**3 . 

T2=V2-V*CREL2 

1 T2 = V*SQRT (I.0“CREL2*«2) . _ _ 

FDX2=V*SREL2/U2 

FDV2=T2/U2 . _„ . _ 

FCV52=TT2**2/U2 

FDXS2=1.0/FNU2~Y*SREL2«*2/U2 

FFV2=FNU2*FDV2*FDVS2 

FPSI2=(T2+V*CX*SREL2/SX)/FNy2 

FMIX2=V*SREL2*FDV2 

_ FTHF2=V*^»FL7^^ ( V2-2^ Q»CRFL2~3- Q*V2*Y* C SRFL2/FNU2 ) **2 )/U2 

CREL3=CXK+SXK*C0S(0.79667*3.1416) 

5RFi_3=SX*CK-“CX*SK*COS CO. 79667*3. 1416) 

FNU3 = S0RT (Vj + Y-2.0*V2*V*CREL3 ) 

FND3=SQRT ( 1.0 + Y“2.0*V«$K«GQS (0^.79667*3.1416) ) 

L3=FNU3**3 

T3 = V2~V*CR£L3 .... _ _^ 

1T3=V*SQRT(i.0-CREL3*«2) 

FDV3=T3/U3 

FDVS3=TT3*«2/U3 - . _ _ . . 

FDXS3=1.0/FNU3-Y*SREL3«*2/U3 

__ FFy 3 =:Fniij3*fl»V3»FDVS3 __ 

FPSI3=(T3+V*CX*SREL3/SX)/FNU3 
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_tL3*FUV3 .. 

FTHE3=V*SRFL3*(V2-2.0*CREL 3^3.C*V2*Y*(SREL3/FNU3)**2)/U3 

__C_REL^=C^K + SXK*CCS(Q^9.dLC2^;»1^1416l . . ... 

SREL4=SX*CK-CX*SK*CUS(C.S6029*3.1416) 

_FNU t = SQRT l.W + Y-2.0*V2*V*CREL4) .. ... . 

FNr)4^SQRT ( 1.0 + Y-2.0*V*SK*C0S (0.96029*3. 1416) ) 

_ L4 = F\ i U4* »3.. _ . . _ ... 

T4 = V2"-V*CRFL4 

TT4 = V*SCRT( I.Q-C RFL4»* 2) . . .. _ _ 

FDX4=V*SREL4/U4 

FDV 4 = T4/U4 _ _ ___ ___ 

FCVS4=rT4**2/U4 

FDXS4=1 .Q/FNU4-Y*SREL^*»2/U4 ^ .. 

FFV4=F^vlU4*FLiV4*FCVS4 

_FPSJ4=(T4 + V*CX*SREL4/SX)/FNU4 __ _ . .. 

FN IX4=V*SREL4*FCV4 

_FTHl4 = V*SRFL4*(V2-2.0*CREL4-3^0*V2*Y*(SRFL4/FJMU.41**2)/.U4.. _ 

GNU-0.36268*FNU1+0.31371*FNU2+0.22238*FNU3+0.10123*FNU4 

__ GND-0 . 36263*FN01+0.3. Lill *Em 2+ 0. 22 23 8»FND3 + 0. 1 Q123» FND4 .. 

GDV-(U 3 626 8*FL:Yi+0.3 1371*FCV2 + 0.22238*FDV3 + 0.10123*FDV4 

G DX^O.36268*FDX1+Q.31371*FDX2+Q.22238*FDX3+Q.10123»FDX4 _ _ 

GDVS-O.3626d*FDVS1+0.31371*FCVS2+0.22236*FDVS3+0.10123*FDVS4 
£0X S-0.36.268*F CXS 1+0.3 13 71*FCXS2 + 0.2223a*FDAS3 + 0.LD123L*fDXS4 
GFV=0.3626R*FFV1+0.31371*FFV2+0.22238*FFV3+0.10123*FFV4 
GP$I -Q. 3A2Jib .*FPS11+0.31371.*EPS 12 + 0.22238 *FPS 13+ 0.10123*FPS14 
GMr<-0.36266*FMlXl+0.31371*FMIX^+0.22238*FMIX3+0.10123*FMIX4 

__GTH..-0.36 268 *F THE 1+0.313 7i*FTHE2+_Q.2223B*ETHE3>0* 10123*F.TJ3E4_ 

NUM l=nuM1 + FLGAT(KTJK(J,K ) )*GNU 

_NUY2=NUY2 + FLnATlKTJK (JiK D *.GDV „ _ ___ . __ 

3-NUM3 + FLCAT( KT JK ( JtK )} ♦GCX 

_N_]^4=NUM.4+FLCAr (KTJK LJiK))*GnVS, . . 

MM^^ = NUM5+FL0AT ( KT JK ( J ,K ) )*GDXS 

_ AUM6-NUM6 + FLQAT(KTJK( J.K ))*GFV .... _ _ , .. 

NUM/=NUM7 + FLCAr(KTJK(J»K ) )*GPSI 

___MJM.3-NUM8 + FLCAT(KJJKt J,K ))*GMIX 

NUM0=MUM9+FLGAT(KTJK(J,K))*GTHE 

_ CENC ML-DENCj’^l + FLCATl KTJK11..X) )*GND . . . 

32 CCNTINUE 

31_CCNfIMUF ... , ....... 

PMAX=RMAXl*\UMl/CENOMl 

_CELTAV(M.L)--4.B0»S»NLM2/NUM1 . . 

CEL r AX ( L )—4.50*S*NLM3/( V2*NUMi) 

_ CELVSQ = 2. 25*S*N UM4/NUM1 __ 

CELTHA=2.25*S*NUM5/(W*NUM1) 

__CVSCAV (M,L ).-SGRT(U_EL^LSlLi„_.. 

CXSuAV(Mt L)-S0RT{0ELTFA) 

_FJ.EPM(M,L)=-6.73*S*NUM6y.NUMl ... „ .. .._ 

CPS I {M,L)=2.25*S*NUM7/(W*V2*SX*SX*NUM1} 

_M.11 =2.^23* S*NUM 8 21 VZ*J1UM ... . „ 

FR rCH( M,L )-2.25*S*NUM9/( Vs*NUMl ) 

_ _ IV^GP tK^LJ.= 0.5* (l*0.+ ( 0ELTAV(M ,L l-FTERMlMtL) +0.3*V2* LDELIHA+SX**2*_ 

ICPSKMtL) ) )/DVSQAV(M,L) ) 

_JLWCLUM^L)=0.5*lJ.^a+ID£LTAXlM,L)-failH(f’,Ll_TaMIX(M*LI/V2L+Q.5*3X*CX* 

ICPST(M,L) )/DXSQAV(M,L) ) 

_ lMAxLM..iJ^IhlXlRMAX+^0..3) .. . . ... _ . , _ , . ... 
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Pfl tflNT INIJF __ 

29 CONTINUE 

Ct]rnR=n^flqAF-ifl*fl MAx*#?»<;QRTf FfiPT) * DENOMl /float rjczoXi__ 

WRITE(6,62A)LL 

hiRITF(fi./i32) _ 

kRITEC6,633)(L,CDVSQAV(M,L),1,11),1,13 ) 

WRTTF(^.fi34) . _ 

WRITE!6,633)(L,!DXSQAV(M,L),M=1,11),L=1,13) 

kR!TEC6e635) _ 

WRITE(6,633)(L,(TW0P(r^,L),M=l,ll),L=l,13) 

kRITFf6,636) ____ 

WRITE!6,6 33)!L,!TWOR!P ,L ),H=1,11),L=1,13) 

WRTTFf6.637) __ 

WRITE!6,638) !L,!PMAX!M,L ), 1, 11) ,L==l, 13 ) 

CO 17 J= 1 e 11 _ 

KPARGV!J)=0 

K J C n =0 __ ___ 

LJ!J)=0 

nni 7K=i, 13_ __ 

KMARGX!K)=0 

17 KTJK(J,K)s=0 _ 

KT0T=0 

nVNSQ=0^0 _ 

SUMNA=0.0 

iBAQ=0____ 

SIB=1,0/SQRT!E0PT) 

r START QF walking AMAx PARTIGLES, ONE AT A TIME ~ _ 

C090A=1,AMAX 

CALL RANDCR) _ 

C NEWTCN RAPHSON DETERMINATION OF V!R) FROM R!V), BUOKER^S DISTRIBUTION 
VC-I-225»SIB«ARSINIR«R} _J__ 

1 2=1.073*V0/SIB 

FV=ERF!Z)-R — _ 

FDV=1.21/!SIB*EXP!Z*Z) ) 

F = FV/Fnv _ __ 

V2=V0-F 

IF lAftS CF) i^.3>2 _____ 

2 V0=V2 

CO TO 1 _ 

3 VINJ=V2 

VHT-0^1 _ 

C TALLY OF INJECTION MARGINAL DISTRIBUTION IN V 

rcAJ^ljll _— 

IFCVINJ.LT.VHT)GO TO 5 

VHT=YHT+Q is 15 _ 

J = ll 

5_ +1 __ 

X2=l*57079633 

K = 1 __ 

2.1 VHH=VHI*TR 

rFLVHH=D£LVHI*TR _-_ 

C LOCATION OF INDICES IN STEP SIZE AND PROBABILITY MATRICES 

MsrFIX!CV2-VHHI/DELVHH+2^0) _—__ 

IFIM.LtIiTgO Tp 18 

IFCMrGTelDGQ TO 20 _—. 

GO TO 19 
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18 


N = 1 _ 

GC ro 19 

20 . _ 

19 L=IFIX( (X2-AHI )/CELXHI + 2.0) 

C .l ANDGM ^ALK Rj?.npEft„ 

IX = 0 

__ IY = J _ _ _ _ 

^^X^M^!AX( M,L ) 

CC2SM,vi^l . „ _ 

CALL RANC(RNR) 

CA LL RANC(R!SjS.) . 

IF( ’N'R.LE.TrtCP(M,L) )GC TC 22 

_ I Y = I Y- 1 __ . 

GC TO 23 
_22. .IY = .IY+l„ 

23 IF(RNS.LE.T^vCR(M,L))GC TL 24 

IX=IX-1 ... _ 

GC TO 25 

24 IX=IX+1 _ 

25 CONTINUE 

C L OCATING TEST PARTICLE IN VELOCITY SPACE 
V2=V2+FL0AT{IY)*DVS0AV(M*L) 

_ CVNS0=DVNSQ+V2*V2-VINJ*VJNJ 

X2 = X2+FLCAT(IX)*DXSOAV(M ,L) 

_ IF ( X2>LT,XK IN>QR,X2«GT.XPAX)GQ 1C, B6. 

VCUT=V2 

__ )(CUT=X2 . . . _ 

XHT^XHI 

_VHT=VHI ... . . _ 

C FIFLC CISTRIBUTICN TALLY 
rC8 IK^1«13 

IF(XOUT.LT.XHT)G0 TO 82 
■■■81 ^_,>HT"XiiT4.D£LXbI . 

K = 13 

82 CCMIIMUE _ . . . „ . 

C08 3J=lt 11 

__ IF ( VCUT>LT>V Hr )GG TO 64. _ 

83 VHT=VHT+DELVhI 

_J=ll . . 

84 KTJK(J,K)=KTJK{J,K)+l 

_IF(N,GT*5C00)GC TO 80 

85 A = N+1 

__ GC ro 2 1__ 

66 SUPNA=SUMNA+FLaAT(N) 

_C^ACCCLNTING FOR FIELD jENERGY CHANGE „ 

TR 1^1.0-0 VN SC* FLOAT! A )/( SUN^N A*FLC AT ( AMAX ) ) 

__ ._IF( iRl-0.0)67^.87,88 

87 TR^O.OOl 

__ GC IG 89_ _ . 

68 TR=SQRT(TRl) 

_.89_ . .. CCNIINUE. ._ . . , . . 

VHT^VHI 

TALLY.CE PARTICLES. EiMTER ING LOSS CC.^IE . 

CC7j=lt11 

_If (.y.2.LT.VhI)aO TO 8 .. 

7 VHT=VHT+OELVhI 
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J = ll _ _ __ _ . . _ 

8 LJ{J)=LJ{J)+l 

GC TO 90 _ — - _ 

80 ieAC=IBAC+l 

90 .. CCNTTNUE 

91 ^K=FL0AT(AMAX-IBAC) 

VSQFVT = 0,0- _ _ _ .. 

C TALLV FOR MARGINAL OISTRieUTION IN V ANO THETA 
C0i6lJ=i9ll - 

CC160K=1,13 

160 KM ARGV C JJ =KHARG\/ ( J ) + KT JK (J •KJ 
V=VHI+(FLCAT(J)-1.5)*CELVHI 

C. CAI Cl.! A r IQN OF. F TFI n FN FRC Y _ _ . .. 

VSGFV(J) = V*V*FLUAT(KMfiRGV(J) ) 

_ VSG.FVT = V-SX|FyT + VSQFVC J } _ . 

161 CONTINUE 

nci63K=1^13 - - _ . „ . 

CG162J=1, 11 

162 KNARGX_(!ll=KMAiLGX.CK) + KTJKiJ-KI __ 

163 CGNIINUE 

CCi64K=i^a3- - - ^ _ 

164 KTC1=KT0T+KHAKGX(K) 

Ti,RITFf 6*4001 . . 

URirE(6t401)(KJ(J),J=l*ll) 

_UR I TF( 6.67^; >. „ . __ 

URITF(6,628)(K,(KTJK(J,K),J=l,ll)t KMARGXIK),K=1,13 ) 
- 4a.R! rEl6^6291 (KMAaGVl J-4-f J=I-tlX) 

V^RITE(6*641 ) (VSQFVIJ ) t J=l, II ) * VSOFVT 

.UR I rEt6.*63QJJ<:TCJ „ .. ... . 

URITF(6,530)QUCCR 

UR ITFC 6.402-1 __ _ . — - -- -- 

UR1TE(6T403)(LJ(J),J=1,11) 

UR.IJE(6.64CJSUKf^m -.. . . _ 

URITEC6,6CC)EOPT♦BMAX,RMAX 1,AMAX,XMIN,XMAX 

30 CCNIINUE . . . _ . 

GC TO 10 

170 S THP „ .. . . _ _ _ 


400 


4 01 
402 


403 

600 

601 

600 


624 

625 


628 

629 


FCRMAT(1HL,46X,41HTALLY CF POINTS IN INJFCTIOil DISTRIBUTION 
!/30-X.3HJ=l .6X * 3HJ^2* 6X * 3hJ-3 * 5X. 3HJ~4.5X-3HJ = 6t 5X.^.3ii-J = 6^ 6X* 

25X,3HJ = 8* 5X*3HJ = 9,4X,4HJ = 10,4X,4HJ = 11) 

FCRMATliQX-5HKJiJ)*10X.illBl _ _ 

FORMAT!1HL*45X,41HTALIY LF POINTS IN LOSS CONE DISTRIBUTION 
1/ 30 X .3-H J= 1 - 6X.*_3HJ=2. 5X - 3F.J=3 - 5X. 3H J = 4 .5X.3HJ=5*1>X- 3N.j = A , sy ^ 3H J = T * 
25X,3HJ=8,5X,3HJ=9,4X,4HJ=10,4X,4HJ=11) 

FOR^ATI lOXs^HLllJJ s-lQX* 1 il 8 ) _ 

FORMAT(F6.1,2!El 4.7),16,13) 

FORMAT (6! FlOeB) ) _ _ 


FORMAT(IHl,60X,23HRAN0nM WALK CALCULATION/IHO,lOX,7HE0•/T =* 

L F6^ 1.5X- 6 HEMAX ^ IP E 10^ 3- 5.X .AHRM.A.Xl^. iPElQ^ 3-3X - 5HAMAX^ , I 5 .6X , 
2 5HXMIN=^0PF10.7,5X,5FXMAX=,FI0*7) 

FCRKA^TtlHL, 65X^2QHlTERAT ION -NUMBER LL~eI14 _ _ 


FORMAT(1HL,42X,48HT0TAL TALLY OF NUMBER OF POINTS IN EACH INTERVAL 

1 /lAX. ihK^15X,3HJ=1^ 5X,3HJ=2,5X^3HJ = 3,5X,3HJ=4,5X,3HJ = 5,5X,3H4=6^ 

2 5X,3HJ=7,5X,3HJ=8,5X,3HJ=9,4X,4HJ=1C,4X,4HJ=11,3X,8HF(THfcTA)) 


FOR Y AT f 1 0 X r I 5 t 10 X ; 12X8) 


FORM AT!llX,4Hf !V),10X, 11 18 ) 
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f.30 FCRVATC 1HL.3X .22H.TGT^L ££.1 NTS-T.ALL I ED . =, J.IO ) _ . . _ . __ 

530 FCRi'.AT(5Xt30HCGLLISION FREQUENCY PAR&METEa=t E14.7) 

C.31 FCR ^AT(6X.4HVhI=.Fm._7^3X^4HXHl^tF10.7,5X.7HD£LVHl = .FlQ.7«5X«7HDEL 

1)IHI = ,F10. 7,5X, 13HPIRRCR RA T 10=, F 11.7 ) 

632 FCR !^AT(1HL.55A,21HSTEP SIZE CVSuAV(F!.L) . . _ . . 

1 /I \ ,2H L. 7X,3hM = 1.8X, 3H^ = 2. dX , 3HM = 3 ♦ 8X , 3H'^ = 4.8X , 3HW=5,3 X , 3HM=6, 

_ 2 8X.3HF^=7.3A» 3 HM = tl,fiX.3HR = 9.7X.4HM-10«7X.4HM = ll)..... 

633 FCRf-ATC IX , 12.1X , 11E1 1 .3) 

634 FCRMAT ( lHL.55XT21hSTEP SIZE CXSQAV(K,L)_.. .. 

1 /Ia,2h L.7X,5HM=1,8X,3HR=2,8X,3HM=3,8X,3hM=4,bX,3HN=5,8X,3HM=6, 

2 8 X.3H V = 7«6X«3FN = 8.3X.3HR=q. 7X«6HiM=10.7X.4HM = il) 

635 FCR.r-'AT( lHL,55X,21HPKaEABILITy TkvOP(K,L) 

_ I /1X.2H L,7X,3FiP=l ,8 X_* 21 8 X t3HMf:.3 t 8X , .3.tiP=Af 8X . 3 HM = 5,8X« 3HM= 6. , 

2 8X,3HM = 7,PX,3HM = 8,8X , >HN = 9,7X14HM=iOt7X,4HM=11) 

636 . FCR"AT(.lHL,55X,21HPR0EA8.ILITY._TW0,k(PiL).. 

1 / 1 ',2H L , 7a, 3HM=1,8X ,3HN‘=2, aX,3HM = 3,8X,3HM=4,6X,3HM = 5,8X,3HM = 6, 

_ 2__8XjJHM=7,HXj.3hM = 8,8X,_3HR = 9,7X_i4HM=lO,7X,4HM = ll). 

637 FCR^ATC lHL,55X,21hFRECUEi\CY «KAX(M,L) 

_ 1 /1A.2H L»7X .3HM=1,8 X , 3HR = 2,8Xj 3_HM.= 3,8X13HM=4j 8X, 3HM=5.j 8 Xj3HM=6.ji, _ 

2 8X,3HH=7,8X,3hP=8,8X,3HP=9,7X.4HP=10,7X,4HM=11) 

FCR^‘AT( lX,.I2,iX.llIll J . .. .... .. 

640 FCRPAK1HL,7FSUP NA=,1 PE 14.7,5X,3HAM = ,1PEI 4.7) 

64J_ _FCRRAT{ llX,6HvVF.(yjt8X,llF8. IjZXjFlO.l) , 

END 



c 


RANDOM WALK CALCULATION 


EQ»/T = 

1.5 fiMAX_^ 1^162F 02 RMAXL= ^.162E OL A«AX= 300 

XMIN« 1.1502463 

XMAX« 1.9913463 

i 

9 

VH= 0.2000000 

XHI* 1.2149463 0ELVHI= 0.2000000 OELXHI= 0.0647000 

MIRROR RATIO* 

1.2GC0168 

L 

8 


_ITERATION _NUHBER LL«1_ 



6 

01 


ti 

Cl 


STFP SIZE OVSOAViMfL) 


L 

f^_=l 

M==2 

M = 3 

M»4 

M=5 

M=6 

K*7__ 

M»8 

H»9 

M-ID 

M*ll 

1 

0.753E-02 

0.7216-02 

C.6466—02 

0.553E-02 

0.459E-02 

0.374E-02 

0.303E-02 

0.245E-02 

0.1996-02 

0.1646-02 

0.1366-02 

2 

0.753E-02 

0.721E-02 

0.646E-02 

0.553E-02 

0.458E-02 

0ja.I4.ErJ12_ 

0.302E-02 

0.245£_r02 

0.199£^02_ 

0.1646-02 

0.1366-02 

3 

0. I58E-02 

0.721E-02 * 

C.646E-C2 

0.552E-02 

0.4586-02 

0.373E-02 

0.302E-02 

0.245E-02 

0.1996-02 

0.1646-02 

0.136E-02 

4 

0.756P-02 

0.721E-02 

C.646E-02 

0.552E-02 

0.458E-02 

0.373E-02 

0.3026-02 

0.2456-02 

0.1996^02 

0.164t-02 

0.136£^02 

5 

0.75HE-02 

0.722E-02 

C .646E—02 

0.552E-02 

0.4586-02 

0.373E-02 

0.302E-02 

0.2456-02 

0.1996-02 

0.1646-02 

0.1376-02 

6 

0.75.3E-02 

0.722E_-02 

0.b46E-02 

0.5526-02 

0.4586-02 

0.373E-02 

Q.3fl?E—0? 0-7456—02 

0.1996-02 

0.1646-02 

0.1376-02 

7 

0.75dE-02 

0.722E-02 

0.6466-02 

0.5526-02 

0.4586-02 

0.373E-02 

0.3026-02 

0.2456-02 

0.1996-02 

0.1646-02 

0.1376-02 

8 

0.75HE-02 

0.7226-02 

0.6466-02 

0.5526-02 

0.458E-02 

0.373E-02 

0. 3026-02 

0.245^^2 

0.1996-02 

0.1646-02 

0.1376-02 

9 

0.75dE-02 

0.722E-02 

0.646E-02 

0.552E-02 

O.450E-O2 

0.373E-02 

0.3026-02 

0.2456-02 

0.1996-02 

0.1646-02 

0.1376-02 

10 

0.75aE~02 

0.7216^02 

0.6466-02 

0.552E-02 

0.4586-02 

0-373E-02 

0.302E-02 

0.245£^02 

0.1996-02 

0.164E-02 

0.1366-02 . 

11 

0.75f3E-02 

0.721E-02 

0.646E-C2 

0.5526-02 

0.4586-02 

0.373E-02 

0.302E-02 

0.2456-02 

0.1996-02 

0.1646-02 

0.136E-02 

_ 12 

0.75dE-02 

0.721E-02 

C.646E-02 

0.553E-02 

0.4586-02 

0.374E-02 

0.3021^02 

0.745E-0? 

0.1996-02 

0.164E-02 

0'.136E-02 

13 

0.7bHE-02 

0.721E-02 

0.6466-02 

0.553E-02 

0.4596-02 

0.3746-02 

0.303E-02 

0.2456-02 

0.1996-02 

0.1646-02 

0.1366-02 


STEP SIZE DXSQAVfH.Ll 


L 

M = 1 

M=2 

M = 3 

H = 4 

M==5 

M*6 

Ms7 

H=8 

M*9 

M*1C 

H«ll 

1 

0.7606^01 

0.2456-01 

C.1388-01 

0.9146-02 

0.6486-02 

Q.48LE-Q2 

Q.369E-Q2 

0.29Q£^Q^ 

Q.233E-Q2 

Q.1916-02 

Q.159E-Q2 

2 

0,7606-01 

0.2456-01 

0.138E-01 

0.913E-02 

0.6486-02 

0.481E-02 

0.368E-02 

0.290E-02 

0.2336-02 

0. 1916-02 

0.159E-02 

3 

0.7596-01 

0.2456-01 

O.l38E-C1 

0.9126-02 

0 . t»47F—02 

Q.4BQE-02 

Q.368E-02 

0.29QE-02 

0.2336-02 

0.1916-02 

0.159E-02 

4 

0.7596-01 

0.2446-01 

C.138E-C1 

0.9126-02 

0.6476-02 

0.480E-02 

0.368E-02 

0.290E-02 

0.2336-02 

0,1916-02 

0.159E-02 

5 

0.7598^01 

0.244E-01 

C.l386-01 

0.912E-02 

0.6476-02 

0.480E-02 

Q.^68E-02 

0.2906-02 

0.233f^Q2_ 

0.I91E-02 

0.159E-02 

6 

0.7596-01 

0.244E-01 

0.1386-01 

0.9126-02 

0.6476-02 

0.4806-02 

0.368E-02 

0.2906-02 

0.2336-02 

0.191E-02 

0.1596-02 

7 

0.759E-01 

0.2446-01 

C.l38E-C1 

0-9116-0? 

0.647E-02 

0.4806-02 

O.360E-O2 

0.290E-02 

0 -2 3 36-0Z 

0.I91F-02 

0.159E-02 

8 

0.7596-01 

0.244E-01 

C.1386-01 

0.9126-02 

0.6476-02 

0.480E-02 

0.368E-02 

0.2906-02 

0,2336-02 

0.1916-02 

0.1596-02 

9 

0.7596-01 

0.244E-01 

C.138E-01 

0.9126-02 

0.6476-02 

0.4806-02 

0.368E-02 

0.2906^02. 

0.2336-02 

0.191E-Q2 

0.1596-02_ 

10 

0.75^6-01 

0.244E-01 

C.1386-01 

0.9126-02 

0.6476-02 

0.480E-02 

0.368E-02 

0.2906-02 

0.2336-02 

0.1916-02 

0.1596-02 

JLL 

0.7596-01 

0.245E-01 

C.l38E-C1 

0-912E-02 

0.6476-02 

0-480E-02 

0.368E-02 

0.2906-02 

0.23 38^0^ 

0.191E-02 

0.159E-02 

12 

0.7606-01 

0.245E-01 

C.1386-01 

0.9136-02 

0.6486-02 

0.481E-02 

0.368E-02 

0.2906-02 

0.2336-02 

0.1916-02 

0.1596-02 

1 3 

0.7b0E-01 

O.245E-0I 

C.1386-01 

0.9146-Q2 

0.6406-02 

Q.4fllE-Q2 

Q.369£^Q2 

0.2906-02 

0.2336-02 

0.1916-Q2 

Q.159E-Q2_ 


PROBABILITY TWOP(M,L) 


Zl 


6 

9 


4 
9 

5 


c 


L 

w = 

1 

Ma. 

2 

M=3 


4 


5 


ft 




n 


1_ 

_M=±C_ 

_H«JJ_ 

1 

0.538E 

00 

0.5116 

00 

C.5C56 

00 

0.5036 

00 

0.5016 

00 

0.500E 

00 

0.5006 

00 

0.4996 

00 

0.4996 

00 

0.499E 

00 

0.4996 

00 

2 

0.5 ^ rf 6 

no 

0.5116 

00 

C.505E 

00 

0.5016 

00 

0.5016 

00 

0.50 IE 

00 

0.5006 00 

0.499£^ 

00 

0.499£_ 

m_ 

_Q^4g9E_ 

00 

0.4996 

00 

3 

0,53-iE 

00 

0.511E 

00 

C.5C5E 

00 

0.5036 

00 

0.5016 

00 

0.501E 

00 

0.5006 

00 

0.499E 

00 

0.4996 

00 

o 

. 

m 

00 

0.4996 

00 

4 

0.5376 

00 

Q.5116 

00 

C.5C5F 

QO 

0.5036 

00 

O.SOLt 

00 

0.50LEl 

00 

0.5006 

00 

Q.499£ 

00 

Q..49.9^ 

m_ 

Q.499A_ 

no 

0.4996. 

no 

5 

0.53 7F 

00 

0.5116 

00 

r.5C56 

00 

0.5036 

00 

0.5016 

00 

0.5016 

00 

0.5006 

00 

0.499E 

00 

0.4996 

00 

0,499£ 

00 

0.499E 

00 

6 

0.5376 00 

0.511E 

00 

0.5056 

00 

0.5036 

09 

0.5016 

00 

0.501E 

00 

0.5006 

00 

0.4996 

00 

0.499£_ 

00 

Q-499E_ 

_QiL_ 

_0.499E_ 

00 _ 

7 

0.5176 

00 

0.511E 

00 

C.5C56 

CO 

0.503E 

00 

0.501E 

00 

0.501E 

00 

0.5006 

00 

0.499E 

00 

0.4996 

00 

0.4996 

00 

0.499E 

00 

8 

0.5 i7F 

00 

0.5116 

00 

C.5C5E 

00 

0.5036 

00 

0.bOLE 

00 

0.501E 

00 

0.500f 

00 

0.499E 

00 

0.499^ 

00 

_Q^4 9 96_ 

00 

0-499^ 

DO 

9 

0.537F 

00 

0.5116 

00 

C.5C56 

00 

0.5036 

00 

0.5016 

00 

0.bOlE 

00 

0.500E 

00 

0.499E 

00 

0.4996 

00 

0.4996 

00 

0.499E 

00 

10 

0.517F 

00 

0.5 1 IE 

00 

C.bC5E 

ilD_ 

0.5036 

J1Q_ 

O.bOLE 

QO 

0.50 IE 

OQ 

Q.5Q0t 

QO 

Q.499£_ 

00 

0 .499^ 

QO 

_Q.499t_ 

QQ 

0.499E 

00 

11 

0.5 1 16 

00 

0.5116 

00 

C.5C5E 

CO 

0.503E 

00 

0.5016 

00 

0.5016 

00 

0.5006 

00 

0.499E 

00 

0.4996 

00 

0.4996 

00 

0.499E 

00 

12 

0.518F 

00 

0.5116 

00 

C.SC5F 

CQ 

0.5036 00. 

0.5016 00 

0.5016 

00 

0.500E 

00 

0.4996 

00 

0.499£_ 

00 

_0.499^ 

_Q0_ 

0.4996 

00 

13 

0.5 136 

00 

0.5116 

00 

C.505e 

00 

0.503E 

00 

0.5016 

00 

0.5006 

00 

0.5006 

00 

0.4996 

00 

0 .4996 

00 

0.4996 

00 

0-.4996 

00 


cn 




PRnnAftM IT Y TUnRfM^I ) 


L 

V = 1 


H=2 


M = 3 


M = 4 


M=5 


M-6 


Ms? 


Ms 8 


Ms9 


M = IC 

Ms^ll 

! 

1 

n - 4 1 ! F 

no 

0 . s n 4 F 

00 

r .. ap?F 

no 

n. '■>01F 

on 

O-sni L 

no 

.n.'iniF _ 

no 

0.^50 OF 

on 

n.snnF 

on 

o^sncF 

■on 

n.cnoF, 

.nn. 

n.5nnF_ 

ilO_ 

9 

2 

O.SQ-JF 

00 

U.5C3C 

00 

C.5C2F 

CO 

0.501E 

00 

0.501E 

on 

0.500E 

00 

0.5ooe 

00 

0.500E 

00 

0.5006 

00 

0.5CCE 

00 

0.5006 

00 

4 

A 

0 . VF 

f)0 

n - ^ F 

on 

r * F 

on 

n , Fo 1F 

on 


no 

O^5i10F 

on 

0. 500F_ 

-Oil- 

n.snnF_ 

on 

n-5onF 

no 

n.cncE- 

.oa. 

0.500E- 

ilO_ 

8 

4 

0.50 •£ 

00 

0.502c 

00 

C.5Ult 

00 

O.SOIE 

00 

0.500E 

00 

0.500E 

00 

0. 500F 

00 

0.500E 

00 

0.5006 

00 

0.500E 

00 

0.500E 

00 

6 


n-4a4F 

00 

n - '=i(M F 

00 

r -sr.i F 

00 

0.ROOF 

00 

a.saoF 

00 

O.SOOfL. 

no 

. O.SOOF no 0.5OOF 

QQ_ 

A-5nOF 

Tia 

Q .FncF 

iKl_ 

n. 5nnF 

on 

01 

6 

C.502P 

00 

n,5Gie 

00 

C,5C0E 

CO 

0.500E 

00 

0.500F 

00 

0.5006 

00 

0.5006 

00 

0.500E 

00 

0.5006 

00 

O.fOCE 

00 

0.500E 

00 

11 

7 

n - «;ntiF 

00 

n.snoF 

on 

r *<^rnF 

on 

._Q^ 500^ 

no 

n.'ifiot. 

on 

_o^ 5noa_ 

no 

n. snoF. 

on 

n.5nnF_. 

■0.0- 

n.5nnF 

JQO_ 

n.cnnF 

.nn 

o-5nnF- 

.00_ 

Zl 

a 

0,49'^F 

00 

0.499t 

00 

C.5C0t 

00 

0.500E 

00 

0.500b 

00 

0.500E 

00 

0.5006 

00 

0.5006 

00 

0.500E 

00 

0.50CE 

00 

0.5006 

00 



n 

nn 

n.4qqF 

on 

P 4 CQF 

on 

n ^ ROOF _ 

on 

n., snnF 

-00 

-fUlOOF on_O^SOOfL 

00 

-Q^SaOE 

■QjQ,- 

fi * p F nn 

p ^ PF 

-on.- 

_0^5nnF_ 

-00_ 


10 

0.49 .F 

ou 

0.498E 

00 

C.499E 

CO 

0.499E 

00 

0.500E 

00 

0.5006 

00 

0.5006 

00 

0.5006 

00 

0.5006 

00 

0.50CE 

00 

0.500E 

00 


11 

n*4q ?<= 

00 

n . 4Q7F 

on 

r *4cqF 

on 

n.499F_ 

00 

n.499t_ 

-Oa 

0.FOO F 

00 

0.5noF 

-00 

o.snnF,, 

lOD- 

n _ 5n PF 

-00 

n^FOCF„ 

-00 

0.5006 

-OQ_ 


12 

0.49 IF 

on 

0.49 ?c 

00 

C.498E 

00 

0.499E 

00 

0.499E 

00 

0.5006 

00 

0.5006 

00 

0.5006 

00 

0.5006 

00 

0.500E 

00 

0.500E 

00 


_L3_ 

_ n^4.-<9K 

-QQ- 

0.496F 

-QQ- 

n -4<;flg_ 

-no- 

0^4QP^ 

on 

0.499.E- 

■■oa... 

0.49.9E- 

-OQ- 

0- 500P. 

■ no 

0.5006- 

jQXL. 

-0...5nn6- 

.nn 

-Q..^5a.CE- 

...nn 

n.500E 

-QQ_ 



FREQUENCY MHAX(M,L) 


1 

^ = 1 


M = '^ 

M.s4 

K1is5 

M=0 

Mss7 

fJiaR 

MsQ 

M= jr 

M»l 1 

1 

22 

23 

25 

21 

30 

33 

37 

41 

45 

49 

54 

y 

yy 


75 

P 2 


_33_ 

_37... 

_4-1_ 

_45 

_50_ 

_54_ 

3 

22 

23 

25 

27 

30 

33 

37 

41 

45 

50 

54 

4 

22 

2 ^ 


? 7 

. Ap 

_3-3 

_ _ 

^ 1 

_45_ 

_ao_ 

_$4_ 

5 

22 

23 

25 

21 

30 

33 

37 

41 

45 

50 

54 

4 

2? 


25 




37 

41 

45 

_so 

54 

7 

22 

23 

25 

27 

30 

33 

37 

41 

45 

50 

54 

Ft 

? ? 

7^ 


2 7 

^p 

3 ' 

37 

41 

45 

5P 

54 - - — 

9 

22 

23 

25 

27 

30 

33 

37 

41 

45 

50 

54 

1 0 

27 

_2"=^ 

25 

27 


3^3- 

37 

_41_ 

_45 

_so 

_54_ 

11 

22 

23 

25 

27 

30 

33 

37 

41 

45 

50 

54 

1 ? 

y y 

2^ 

2*^ 

2 7 

■^n 

^ 4 

^7 

41 


5n 

5^ 

13 

22 

23 

25 

27 

30 

33 

37 

41 

45 

49 

54 


-TALLY OF POINTS IN t-NJECTION Dl-STR IBUT1 ON- 

J=1 J = 2 J = 3 J=A J~5 J=6 J=7 J=8 J=9 3=^10 J^U 

JLJUJ_22_92_9J_ (x5 _ 6Q _^6_10_19-15_2_6. 


K 

J = l 

Js2 

TOTAL TALLY OF 

_la^ ls4 

NUMBER OF 
Js5 

POINTS 

IN FACH 

—4=7 

INTERVAL 

J=A . 


J=10 

J*T 1 

FITHFTAT 

1 

17 

61 

143 

184 

224 

191 

265 

118 

7 

0 

37 

1247 

? 

54. 

1 1 1 

- 2ftfl 

42 3 

525 


7n2 

?74 

- - - 53 

4 

IfQ 

3 1 P4 

3 

25 

200 

434 

745 

889 

989 

782 

428 

76 

104 

755 

4927 

4 

?a 

251 

4qA 

1 n55 

1 1 25 _ 

1 5P4 

772 

431 

1 71 

_ 272 


4450 

5 

£8 

269 

650 

1204 

1215 

1900 

1462 

486 

182 

323 

589 

B366 

A 

42 

2q7 

_ A2fl 

14iq 

1 54] 

2337 

1 704 

flP7 

2fin 

349 

4n9 

inn55 

7 

127 

365 

1041 

1447 

1960 

2352 

1859 

1249 

4S3 

650 

230 

11773 

fl 

4iA 

^26 


) 52 ^ 

1 fl57 

1 

] 3Rn 

1 in<J 

5#R 

49 3 

1^1 

_ 10009 

9 

91 

324 

635 

1234 

1505 

1532 

1044 

704 

6 39 

328 

162 

8198 

-. in 


2y^ 

4q2 

91 1 

) 1 


4fl4 

529 

507 

232 

1 

4442 

11 

21 

168 

313 

667 

866 

1224 

610 

427 

2 99 

299 

99 

4993 

_L2_ 

^9 

11*^ 

247 

3RR 

5 0 5 

8*^0 

- 4f12 

246— - 

1 €3 

200 

36 

040 j. 

‘ 13 

12 

68 

133 

177 

252 

387 

2 75 

104 

99 

U2 

65 

1684 

‘=TV) 

64^ 

277® 

A5^7 

1 132^7. 

_1 3741 

17pn2_ 

1 2P25 

_AQ1.2_ 

35_C7 

3304 

2471 


VVF(V) 

6.4 

250.1 

1634.2 

5574.7 

11146.4 20572.4 

20322.2 

15552^0 10395.3 

12223.5 

11779.1 

109,456.5 


■XCTM.. POINTS TALLIED =_ aQ69Q_ 

CCLLISION FREQUENCY PARAMETER^ 0.1454641E-12 




TALLY 

QF POINTS 

IN LOSS 

CONE 

DISTRIBUTION _ 




» 


Jsl J=2 

4. 

II 

UJ 

J*4 

J = 5 

J=6 

J=7 Js8 

J»9 

J-10 

J*ll 

t 

LJtJ) 

_ 128, _L2J_ 

_II_ 

69 

_ 

33 

_13_9_ 

_4__- 

_1_ 

_1_ 

t 

B 

SUM NA* 7.61880COE 04 

AP» 4.9899999E 02 









§ 



APPENDIX F 


ANALYTICAL SOLUTION FOR CASE OF SHORT WALKS 

Consider the walks so short that the Fokker-Planck coefficients are constant over 
the distance traveled. Assume that the walk terminates when a particle first reaches 
a prescribed 6 distance from its initial location. Effects of V are thus of second 
order and can be neglected. For a spherical field distribution, both equations (4) 
and (5) in combination with a source term as in equation (6) reduce to 

1 <(A0)^> — 4in 0 = -s sin 9 (FI) 

2 d9 \ 39/ 

If the initial test-particle location is at 6 ^ 7r/2, then sin @ « 1 and s = - (7r/2)] 

reducing equation (FI) to 


1 ((A0)2> i?l = -s 6 fe - (F2) 

2 9^2 " I 2 ) 

This model may simulate the end loss problem for mirror ratios very close to 1. 0. 
Injection would be normal to the B field at a constant rate s The boundary conditions 
are 


f(0^) = f(7r - = 0 


(F3) 


Using (FI) in (F2) and integrating both sides of the equation gives 


39 39 



where H is the Heaviside unit function. Integrating a second time and using equation 
(F3) yields 
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m - ^ 

30 


(0 - 0 ^) = 0 


0=0 


2s 


((A0)^> 


2-/0 

.^2^ V 2/ 


when 0 < 0 < - 

^ 2 


when —< 9 ^ (ir - 9 ) 
2 


r (F4) 
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At 9-71-9 


30 


o 


0=0p <(A0)^> 


(F5) 


Using equation (F5) in (F4) gives 


f(0) =-2_ (0 - 0 ) 

<(A0)2> 


<(A0)2> 


when 0 „ ^ 0 ^ — 
^ 2 






► 

—— (tt - 0 - 9) when - < 9 < n <9 
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(F6) 


Loss rate must equal injection rate for steady state so 


/ 7r-0 


s(0)sin 0 d0 = s 



77-9 


9 


(' 


- 


6(0-jsin 0 d0 = s 


2j 


The number density, distribution function relation must be 


% ■ 



77-9 


f(0)sin 0 d0 = 


2s 


<(A0)^) 


0 



77/2 


(0 - 0^)sin 0 d0 


0 


2s 


<(A0)2> 


-— (1 - sin 0^) 


(F7) 
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Using equation (F6) gives 


n = 


nb((A0)^> 

2(1 - sin e^) 


For small - - 9 


c> 


sin 9 


...(Lli 


(F8) 


to second order so that the loss rate becomes 


^ n^^<(A0)^) 

n =- 


(F9) 


- - 9 


Evaluating ((A0) ) for a spherical Maxwellian field distribution (eq. (11c)) and sub¬ 


stituting into equation (F9) yields for v = v 


o 








In /3 

1 

(2kT^)^/2 









ertM-2 


(FIO) 

If the particles are, for example, injected at 10 times the average field energy then 


1 mv2 = M „ 

2 “ 2 


or 


E 


= 15 
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and for deuterons, for example, 



n^log/3 


0.773x10'^® 


(Fll) 


This result is used in figure 6. 

It is interesting to compare equation (F9) with equation (33) of reference 8. The 
result of reference 8, derived strictly from a random walk, is for only one absorbir^ 
wall. The loss rate predicted by equation (F9) (for two absorbing walls) is just twice 
that of reference 7. 

To determine the marginal distribution in 9 for use on figure 7(a) substitute equa¬ 
tions (F7) and (F8) into (F6). This yields 


F(,(9) 


{.( 9 ) e - 






when 9 „ — Q — — 

2 


y ( 17 ) 


when - (9 < (it - 9„) 

2 ^ 
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APPENDIX G 


SELECTION OF RANDOM NUMBERS FROM THE NONUNIFORM 
INJECTION DISTRIBUTION OF REFERENCE 5 

The simulation procedure of sampling from a nonuniform distribution when the com¬ 
puter library contains only a uniformly distributed set of random numbers is quite com¬ 
mon (refs. 9, 21, and pp. 252-264 of ref. 15). The relation between a uniformly dis¬ 
tributed random number R and an arbitrarily distributed random number v is 

/ R 

dr = / f(v)dv 

U •'O 

where it is necessary only that the second integral be a monotone increasing function of 
if. This says that the cumulative distribution function R of the probability density f(v) 
is uniformly distributed in Y. If f(v) is integrable, can be found. But to find 

'f{R) explicitly in the case of interest herein required a root finding method such as, for 
example, the Newton-Raphson method (ref. 22). 

The injection distribution of reference 5 is 



^ o -mvV2kT, 9 
((A^)^) e ^ dv 


R{r) = - — 



9 -mv'^/2kTy. 9 

((A^)^>e ^ dv 


where ((A0)^) is given by equation (11c). Letting x = \/^/2kT^ v results in a new 

expression for 
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ym/2kT^ Y 


R(^ = 


0 



OO 


0 




r-Le-\( 

X - — jerf X 


^ 2x/ J 

J2 / 1 

\ 1 2 

+ fx-± 

V 2x 

jerf X e”^ c 


-x^ dx 


(Gl) 


This was integrated by the method of Gaussian quadratures. The curve fit 


R(f') = erf (1.073 


VS 


is a good approximation to the result. This is shown on figure 12 where V = 



Figure 12. - Random number distribution to fit injection distribution of 
reference 5. 
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I 


used as the abscissa and \/m 72 kT^ was set equal to V3/2 to make ref¬ 

erence velocity equal to the mean field velocity. 

For an initial approximation V(R) for use in the Newton-Raphson method, the curve 


V = arcsin R 


2 


was used. 

Results of this procedure to generate the injection distribution is shown to be satis¬ 
factory in figure 8(c). 
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